{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAO0AAABECAYAAABtXrKpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAB/klEQVR4nO3YMYrUYBzG4X+GDaLDuIWewGIu4BksvIqVILZewMrGUvAUFp7BC0xrJ6uwu2RFAvNZyFo4RhA2m331ecp8zRvIj5mka621AmKslh4A/B3RQhjRQhjRQpijqYP9fl/DMFTf99V13XVugv9aa63Gcaz1el2r1eHv6mS0wzDUbrebdRwwbbvd1mazObg+GW3f91VV9ezVmzo5PZtv2YLevnhar989X3rGbJ48flmr9x+WnjGb/aOHdbF/sPSMK9fVWLdXH382+KvJaC//Ep+cntWnL6fzrLsBTi8+Lz1hVt3Xb0tPmFWr3z/Y/4Kp11IfoiCMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCGMaCHM0dRBa62qqu4f3722MUs4vnNv6QmzardvLT1hVl2NS0+4cpf3dNngwXmbODk/P6/dbjffMuCPttttbTabg+uT0e73+xqGofq+r67rZh8I/NBaq3Eca71e12p1+AY7GS1wM/kQBWFEC2FEC2FEC2G+A7O0WWjJZuMgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 288x72 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# load required packages\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns; sns.set()\n",
    "import statsmodels.api as sm\n",
    "import pickle\n",
    "import joblib\n",
    "from __future__ import division\n",
    "from IPython.display import display\n",
    "\n",
    "# set styling for plots\n",
    "%matplotlib inline  \n",
    "sns.set(style = \"whitegrid\")\n",
    "sns.set_palette('cubehelix',4)\n",
    "plt.rc('text', usetex=True)\n",
    "plt.rc('font', family='serif')\n",
    "sns.palplot(sns.color_palette())\n",
    "\n",
    "# define names of models\n",
    "clfs = [\"LogitNonLinear\",\"RandomForestIsotonic\"]\n",
    "outnames = {\"LogitNonLinear\": \"Nonlinear Logit\",\"RandomForestIsotonic\" : \"Random Forest\",\"xgboost_output\": \"XGBoost\",\"Difference\" : \"Difference\"}\n",
    "\n",
    "# set race categories for plots\n",
    "plotrace = [\"Asian\",\"White Non-Hispanic\",\"White Hispanic\",\"Black\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load and clean data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set paths to saved data\n",
    "path = '../../data/' \n",
    "\n",
    "# load main dataset \n",
    "allvals = pd.read_csv(path + 'all_vals_race1_interestrate1.csv',index_col=0)\n",
    "\n",
    "# load predictions \n",
    "preds00 = pd.read_csv(path + '../output/_race0_interestrate0.csv',index_col=0) \n",
    "\n",
    "# load saved classifiers and save feature names\n",
    "clfpath = '../../output/'\n",
    "models = {}; features = {}\n",
    "for name in clfs:\n",
    "    if name.startswith('Logit'):\n",
    "        models[name] = pickle.load(open(clfpath + name + \"_race0_interestrate0.pkl\",'rb'))\n",
    "        features[name] = list(models[name].params.index)\n",
    "    else:\n",
    "        models[name] = joblib.load(open(clfpath + name + \"_race0_interestrate0.pkl\",'rb'))\n",
    "        features[name] = [x.lower() for x in list(pd.read_csv(path + 'feature_names_norace' + name + '.csv').columns)]\n",
    "features['RandomForestIsotonic'].remove('sato')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# clean race variable in main data\n",
    "racecols = [col for col in list(allvals) if col.startswith('race_dum')]\n",
    "allvals[\"White Non-Hispanic\"] = 1 -  allvals[racecols].sum(axis=1)\n",
    "allvals[\"Race\"] = allvals[racecols + [\"White Non-Hispanic\"]].idxmax(axis=1).str.replace('race_dum_','').replace('White hisp','White Hispanic')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# merge main data and predictions \n",
    "#df0 = allvals.drop(['LogitNonLinear','RandomForest','RandomForestIsotonic'], axis=\"columns\")\n",
    "df0 = allvals\n",
    "df0 = df0.join(preds00)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Figure 4: Example heatmaps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# utility for adding fico, ltv and income bins used by logit to a dataframe\n",
    "def mkbins(df):\n",
    "    bins = pd.DataFrame(index = df.index)\n",
    "\n",
    "    fico_cuts = [0] + list(range(280,870,20))\n",
    "    fico_bin = pd.cut(df[\"fico_orig_fill\"], fico_cuts, labels = fico_cuts[0:-1], right = False).fillna(0)\n",
    "    fico_bin.loc[(fico_bin>0) & (fico_bin < 600),] = 600\n",
    "    fico_bin.loc[fico_bin==840,] = 820\n",
    "\n",
    "    ltv_cuts = list(range(20,110,5))\n",
    "    ltv_bin = pd.cut(df[\"ltv_ratio_fill\"], ltv_cuts, labels = ltv_cuts[0:-1], right = False)\n",
    "    ltv_80 = (df[\"ltv_ratio_fill\"]==80.0).astype(int)    \n",
    "\n",
    "    inc_cuts = list(range(-25,550,25))\n",
    "    income_bin = pd.cut(df[\"applicant_income\"],inc_cuts, labels = inc_cuts[0:-1],right = False)\n",
    "\n",
    "    bins = bins.join(pd.get_dummies(fico_bin, prefix = \"fico_bin_dum\"))\n",
    "    bins = bins.join(pd.get_dummies(income_bin, prefix = \"income_bin_dum\"))\n",
    "    bins = bins.join(pd.get_dummies(ltv_bin, prefix = \"ltv_bin_dum\"))\n",
    "    bins[\"ltv_80_dum_0\"] = 1 - ltv_80\n",
    "    bins['const'] = 1\n",
    "    return bins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=2)]: Using backend ThreadingBackend with 2 concurrent workers.\n",
      "[Parallel(n_jobs=2)]: Done   2 out of   2 | elapsed:    0.2s remaining:    0.0s\n",
      "[Parallel(n_jobs=2)]: Done   2 out of   2 | elapsed:    0.2s finished\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-13-ec8765fd29a4>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     50\u001b[0m \u001b[0mheat\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"FICO\"\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcut\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mallvals\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"fico_orig_fill\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mfic_grid1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mlabels\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfic_grid1\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'int'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     51\u001b[0m \u001b[0;31m# heat[\"Minority\"] = ((race == \"Black\") | (race == \"White Hispanic\")).astype('int')\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 52\u001b[0;31m \u001b[0mheat\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"Minority\"\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mrace\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m\"Black\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'int'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     53\u001b[0m \u001b[0mheat\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"Majority\"\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mrace\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m\"White Non-Hispanic\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'int'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     54\u001b[0m \u001b[0minrange\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mheat\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"Income\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnotnull\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m&\u001b[0m \u001b[0mheat\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"FICO\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnotnull\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/anaconda3/lib/python3.8/site-packages/pandas/core/frame.py\u001b[0m in \u001b[0;36m__setitem__\u001b[0;34m(self, key, value)\u001b[0m\n\u001b[1;32m   3161\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3162\u001b[0m             \u001b[0;31m# set column\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 3163\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_set_item\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   3164\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3165\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_setitem_slice\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkey\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mslice\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/anaconda3/lib/python3.8/site-packages/pandas/core/frame.py\u001b[0m in \u001b[0;36m_set_item\u001b[0;34m(self, key, value)\u001b[0m\n\u001b[1;32m   3240\u001b[0m         \"\"\"\n\u001b[1;32m   3241\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_ensure_valid_index\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvalue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 3242\u001b[0;31m         \u001b[0mvalue\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_sanitize_column\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   3243\u001b[0m         \u001b[0mNDFrame\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_set_item\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkey\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3244\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/anaconda3/lib/python3.8/site-packages/pandas/core/frame.py\u001b[0m in \u001b[0;36m_sanitize_column\u001b[0;34m(self, key, value, broadcast)\u001b[0m\n\u001b[1;32m   3874\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3875\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvalue\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mSeries\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 3876\u001b[0;31m             \u001b[0mvalue\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mreindexer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvalue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   3877\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3878\u001b[0m         \u001b[0;32melif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvalue\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mDataFrame\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/anaconda3/lib/python3.8/site-packages/pandas/core/frame.py\u001b[0m in \u001b[0;36mreindexer\u001b[0;34m(value)\u001b[0m\n\u001b[1;32m   3860\u001b[0m                 \u001b[0;31m# GH 4107\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3861\u001b[0m                 \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 3862\u001b[0;31m                     \u001b[0mvalue\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreindex\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mindex\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_values\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   3863\u001b[0m                 \u001b[0;32mexcept\u001b[0m \u001b[0mValueError\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0merr\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3864\u001b[0m                     \u001b[0;31m# raised in MultiIndex.from_tuples, see test_insert_error_msmgs\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/anaconda3/lib/python3.8/site-packages/pandas/core/series.py\u001b[0m in \u001b[0;36mreindex\u001b[0;34m(self, index, **kwargs)\u001b[0m\n\u001b[1;32m   4343\u001b[0m     )\n\u001b[1;32m   4344\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mreindex\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mindex\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 4345\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0msuper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreindex\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mindex\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mindex\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   4346\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4347\u001b[0m     def drop(\n",
      "\u001b[0;32m~/opt/anaconda3/lib/python3.8/site-packages/pandas/core/generic.py\u001b[0m in \u001b[0;36mreindex\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m   4809\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4810\u001b[0m         \u001b[0;31m# perform the reindex on the axes\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 4811\u001b[0;31m         return self._reindex_axes(\n\u001b[0m\u001b[1;32m   4812\u001b[0m             \u001b[0maxes\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlevel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlimit\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtolerance\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfill_value\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcopy\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4813\u001b[0m         ).__finalize__(self, method=\"reindex\")\n",
      "\u001b[0;32m~/opt/anaconda3/lib/python3.8/site-packages/pandas/core/generic.py\u001b[0m in \u001b[0;36m_reindex_axes\u001b[0;34m(self, axes, level, limit, tolerance, method, fill_value, copy)\u001b[0m\n\u001b[1;32m   4825\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4826\u001b[0m             \u001b[0max\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_get_axis\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 4827\u001b[0;31m             new_index, indexer = ax.reindex(\n\u001b[0m\u001b[1;32m   4828\u001b[0m                 \u001b[0mlabels\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlevel\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mlevel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlimit\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mlimit\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtolerance\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtolerance\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmethod\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4829\u001b[0m             )\n",
      "\u001b[0;32m~/opt/anaconda3/lib/python3.8/site-packages/pandas/core/indexes/base.py\u001b[0m in \u001b[0;36mreindex\u001b[0;34m(self, target, method, level, limit, tolerance)\u001b[0m\n\u001b[1;32m   3523\u001b[0m                             \u001b[0;34m\"with a method or limit\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3524\u001b[0m                         )\n\u001b[0;32m-> 3525\u001b[0;31m                     \u001b[0mindexer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmissing\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_indexer_non_unique\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtarget\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   3526\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3527\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mpreserve_names\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mtarget\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnlevels\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m1\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mtarget\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mname\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/anaconda3/lib/python3.8/site-packages/pandas/core/indexes/base.py\u001b[0m in \u001b[0;36mget_indexer_non_unique\u001b[0;34m(self, target)\u001b[0m\n\u001b[1;32m   4941\u001b[0m             \u001b[0mtgt_values\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtarget\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_get_engine_target\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4942\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 4943\u001b[0;31m         \u001b[0mindexer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmissing\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_engine\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_indexer_non_unique\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtgt_values\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   4944\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mensure_platform_int\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mindexer\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmissing\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4945\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32mpandas/_libs/index.pyx\u001b[0m in \u001b[0;36mpandas._libs.index.IndexEngine.get_indexer_non_unique\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32m<__array_function__ internals>\u001b[0m in \u001b[0;36mresize\u001b[0;34m(*args, **kwargs)\u001b[0m\n",
      "\u001b[0;32m~/opt/anaconda3/lib/python3.8/site-packages/numpy/core/fromnumeric.py\u001b[0m in \u001b[0;36mresize\u001b[0;34m(a, new_shape)\u001b[0m\n\u001b[1;32m   1428\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1429\u001b[0m     \u001b[0mrepeats\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mnew_size\u001b[0m \u001b[0;34m//\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msize\u001b[0m\u001b[0;34m)\u001b[0m  \u001b[0;31m# ceil division\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1430\u001b[0;31m     \u001b[0ma\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mconcatenate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mrepeats\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mnew_size\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1431\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1432\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mreshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnew_shape\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<__array_function__ internals>\u001b[0m in \u001b[0;36mconcatenate\u001b[0;34m(*args, **kwargs)\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "# make predictions on fico/income grid in figure\n",
    "\n",
    "# grid of income / fico\n",
    "inc_grid = np.linspace(20,200,50)\n",
    "fic_grid = np.linspace(650,820,50)\n",
    "inc_list = [x for x in inc_grid for y in fic_grid]\n",
    "fic_list = [y for x in inc_grid for y in fic_grid]\n",
    "\n",
    "# set up grid\n",
    "G = len(inc_grid)*len(fic_grid)\n",
    "K = len(features[\"RandomForestIsotonic\"])\n",
    "x = pd.DataFrame(np.zeros((G,K)),columns = features[\"RandomForestIsotonic\"])\n",
    "\n",
    "# insert income/fico grids\n",
    "x[\"applicant_income\"]=inc_list\n",
    "x[\"fico_orig_fill\"]=fic_list\n",
    "\n",
    "# hard code the remaining variables as explained in paper \n",
    "ltv = 79.99 # exactly 80 doesn't work\n",
    "amt = 300000\n",
    "x[\"ltv_ratio_fill\"] = ltv\n",
    "x[\"orig_amt\"]=amt \n",
    "x[\"log_orig_amt\"]=np.log(amt)\n",
    "# x[\"sato\"] = 0\n",
    "x[\"occupancy_type_dum_1\"]= 1\n",
    "x[\"investor_type_dum_2\"] = 1\n",
    "x[\"loan_purpose_dum_1\"] = 1\n",
    "x[\"orig_year_dum_2011\"] = 1\n",
    "x[\"prop_state_dum_CA\"] = 1\n",
    "x[\"prop_state_dum_TX\"] = 0\n",
    "x[\"cur_int_rate\"] = 4.5\n",
    "\n",
    "xbins = mkbins(x)\n",
    "x = x.join(xbins)\n",
    "\n",
    "# predictions for PD contour plot \n",
    "pred = pd.DataFrame(index = x.index, columns = clfs)\n",
    "for name in clfs:\n",
    "    if name.startswith('Logit'): \n",
    "        pred[name] = models[name].predict(x[features[name]]) \n",
    "    else:\n",
    "        pred[name] = models[name].predict_proba(x[features[name]])[:,1]\n",
    "        \n",
    "# get group densities on grid\n",
    "race = df0['Race']\n",
    "inc_grid1 = np.linspace(20,200,10)\n",
    "fic_grid1 = np.linspace(650,820,10)\n",
    "heat = pd.DataFrame(index = allvals.index, columns = [\"Income\",\"FICO\",\"Minority\"])\n",
    "heat[\"Income\"] = pd.cut(allvals[\"applicant_income\"],inc_grid1,labels = inc_grid1[1:].astype('int'))\n",
    "heat[\"FICO\"] = pd.cut(allvals[\"fico_orig_fill\"],fic_grid1,labels = fic_grid1[1:].astype('int'))\n",
    "# heat[\"Minority\"] = ((race == \"Black\") | (race == \"White Hispanic\")).astype('int')\n",
    "heat[\"Minority\"] = (race == \"Black\").astype('int')\n",
    "heat[\"Majority\"] = (race == \"White Non-Hispanic\").astype('int')\n",
    "inrange = heat[\"Income\"].notnull() & heat[\"FICO\"].notnull()\n",
    "heat = heat[inrange]\n",
    "piv = {}\n",
    "piv['Black'] = pd.pivot_table(data = heat, index = \"FICO\", columns = \"Income\"\n",
    "                            , values = \"Minority\", fill_value = 0,\n",
    "                            aggfunc = (lambda x: x.sum()/heat[\"Minority\"].sum()))\n",
    "piv['White Non-Hispanic'] = pd.pivot_table(data = heat, index = \"FICO\", columns = \"Income\"\n",
    "                            , values = \"Majority\", fill_value = 0,\n",
    "                            aggfunc = (lambda x: x.sum()/heat[\"Majority\"].sum()))\n",
    "\n",
    "# draw contour plots\n",
    "X,Y = np.meshgrid(inc_grid,fic_grid)\n",
    "fig,ax = plt.subplots(2,len(clfs),figsize = (10,8))\n",
    "ext = [inc_grid1.min(),inc_grid1.max(), \n",
    "       fic_grid1.min(),fic_grid1.max()]\n",
    "cvals = [[0.1,0.2,0.3,0.5],[0.2,0.3,0.6,1,1.3,1.6]]\n",
    "for (i,name) in enumerate(clfs):\n",
    "    for (j,g) in enumerate(piv.keys()[::-1]):\n",
    "        # dist for group j\n",
    "        ax[i,j].imshow(piv[g][::-1],cmap=\"Purples\",extent = ext,alpha = 0.6)\n",
    "        \n",
    "        # contour for clf i overlaid\n",
    "        Z = pred[name].values.reshape(X.shape)\n",
    "       \n",
    "        CS = ax[i,j].contour(X,Y,100*Z,cvals[i],cmap=None,colors='k')\n",
    "        plt.clabel(CS, fontsize=14,inline=True,fmt='%1.1f')\n",
    "        \n",
    "#         ax[i,0].set_ylabel('%s' % (outnames[name]),fontsize=18)\n",
    "        ax[i,j].set_xlabel('Income',fontsize = 14)\n",
    "        ax[i,j].set_ylabel('FICO',fontsize = 14)\n",
    "\n",
    "\n",
    "        ax[i,j].set_title('%s: %s' % (outnames[name],g),fontsize=15)\n",
    "\n",
    "\n",
    "        ax[i,j].grid(False)\n",
    "        ax[i,j].tick_params(labelsize = 16)\n",
    "\n",
    "\n",
    "plt.tight_layout(pad=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Figure 5: CDFs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAewAAAFPCAYAAACPov0gAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA2JElEQVR4nO3dfWwcZ34f8K+cs+SztbM8pWibmkMU8cEXa5d3TWX2yhGaApYgzjIpUtGShgVygWhzqcCARQRYOgki0ueVgCTiBgkVIDhyJEvIBQhHzvmaNtAuhThoA+zQsY59CXfIxogTgEMhTdpanFmeY/lt+wc94519n+W+kt8PIIg7Mzvz29md/e3zzPNyIJfL5UBEREQd7ZF2B0BERETVMWETERF1ASZsIiKiLsCETURE1AWYsImIiLoAEzYREVEX+FKrD2jbNjRNw9bWFqampqpur6oqRFGEZVkAAEVRmh0iERFRx2lpCVvXdei6DtM0kc1mq24/OzsLURQhyzIURcHGxgZSqVQLIiUiIuosLU3YkiRBlmUEAoGatr99+zZkWXYfDw8PQ9O0ZoVHRETUsTr2HrZhGEXLBEGAruttiIaIiKi9Wn4Pu1aWZSEYDHqWFT6uxWeffYYf/vCHePTRR3HgwIFGhUdERFSTXC6Hjz/+GE888QQeeaT+cnLHJmzbtiuuEwShpv388Ic/xLvvvtuosIiIiOry9NNP13xLuJSOTdiCILgtwx2Fj2vx6KOPAtg5UQcPHmxIbNRYmUwG4XC43WFQGS+99BK2t7fxe7/3e3U9f+k//Rd85Z//K4SeyuGJJ55ocHQE8BrqdB999BHeffddNx/Vq2MTdjAYLCplO49rLV0DcKvBDx48iEOHDjUuQGoovjedy7IsZLPZut+jTz/5DDk8ikcfzfF9biKe286329uyHdvoLBQKFSVmy7IgSVKbIiIiImqfjkrYpml6+llHIhHP43Q6zYFTiIhoX2ppwjYMA6qqYmlpCbquQ1VVT/etVCrl6Wcdj8dhmiZ0XYemaejr6/P0yyYiItovWnoPOxQKIRQKIRqNllwfjUaL1pXbloiIaD/pqCpxIiIiKo0Jm4iIqAswYRMREXUBJmwiIqIuwIRNRETUBZiwiYiIugATNhERURdgwiYiIuoCTNhERERdgAmbiIioCzBhExERdQEmbCIioi7AhE1ERNQFmLCJiIi6ABM2ERFRF2DCJiIi6gJM2ERERF2ACZuIiKgLMGETERF1ASZsIiKiLsCETURE1AWYsImIiLoAEzYREVEXYMImIiLqAkzYREREXYAJm4iIqAswYRMREXUBJmwiIqIuwIRNRETUBZiwiYiIugATNhERURdgwiYiIuoCTNhERERdgAmbiIioCzBhExERdQEmbCIioi7AhE1ERNQFmLCJiIi6ABM2ERFRF2DCJiIi6gJM2ERERF2ACZuIiKgLMGETERF1gS+146CqqkIURViWBQBQFKXq9oIgAABs20Y0Gm16jERERJ2k5SXs2dlZiKIIWZahKAo2NjaQSqXKbq+qKqLRKBRFgaIokCQJqqq2MGIiIqL2a3nCvn37NmRZdh8PDw9D07Sy2yeTSc/jUCiE1dXVpsVHRETUiVqasA3DKFomCAJ0XS/7nGAwiIsXL7qPNU3D8PBwU+IjIiLqVC1N2JZlIRgMepYVPi4Uj8exvLyMgYEBqKqKYDDoKaETERHtBy1tdGbbdsV1TsOyfKIoYmJiArquI5FIYHx8vK6EnclkfD+HWmdlZaXdIVAZ2WwWQP3v0aeffAIA+J9/8Rf48mOPNSwu8uI1tPe1NGELguC2DHcUPi40MzMDRVEQjUah6zomJydhmiauXbvm69jhcBiHDh3yHTM138rKCo4dO9buMKiMQCCAbDZb93u08e7fAwC+8fWvIxAINDI0+hyvoc728OHDhhQaW1olHgwGi0rZzuNSpWvDMBAIBBAKhQAAkiThrbfewvLycvODJSIi6iAtTdihUKgoMVuWBUmSSm5vWRZ6eno8ywRBwODgYLNCJCIi6kgt79YViUQ8/a7T6bRn4BTTNN31kiQVtSC3bRuiKLYmWCIiog7R8pHO4vE4VFWFruswTRN9fX2eRmSpVAq6rrvL4vE4Zmdn0dfX524zNTXV6rCJiIjaqi1Dk1YaWjQajXrWi6LIBE1ERPseJ/8gIiLqAkzYREREXYAJm4iIqAswYRMREXUBJmwiIqIuwIRNRETUBZiwiYiIugATNhERURdgwiYiIuoCTNhERERdgAmbiIioCzBhExERdYG6EvbS0hJeeOEFnDp1yl32xhtvNCwoIiIi8vI9W9ft27exsbGBWCwG0zTd5WfPnsXdu3c9SZyIiIgaw3fCDgaDiMViAADLsjzrcrlcY6IiIiIiD99V4oIglF1n2/augiEiIqLSfCdswzCwvr5etHx9fb2oxE1ERESN4btKfHx8HJOTk1hfX0cgEIAoijBNEz09Pbhx40YzYiQiItr3fCdsAJibm4NpmlheXoZt21AUBYODg42OjYiIiD7nO2FPTk7i61//Ol588UWIotiMmIiIiKiA74R9/PhxRCKRkuu2t7dx+PDhXQdFREREXr4bnYmiWLZxmaZpuw6IiIiIivkuYSeTSRiGgWw2C1EUEQgEAOx06TIMAy+++GLDgyQiItrvfCfsTCaDqakpBIPBonXz8/MNCYqIiIi8fCfsK1eu4OjRo55l2WwWW1tbuHDhQsMCIyIioi/4voddmKwBIBAIoKenB5ubmw0JioiIiLzq6oe9ubnp9sHOt7Gxwck/iIiImsB3wl5eXoaqqjh69Chs23bHFrdtG6+88krDAyQiIqI6Evba2hpef/11AHCn13QGUFlfX8czzzzTwPCIiIgIqOMedm9vr/u3KIp4++233cdbW1sNCYqIiIi8fCdsx927dwEAq6ur2N7eBrBT+iYiIqLG852wh4aGcP36ddy5cwcAEI1G8dxzz+Gb3/wmS9hERERNUlcr8fHxcfdvURTxzjvvYG1trWSXLyIiItq9uqvECx09ehTLy8uN2h0RERHlqauEvby8XDRIim3buHPnDr73ve81JDAiIiL6gu+EnUgkYJpmybmws9lsQ4IiIiIiL98Ju7+/H7FYrOw6IiIiajzf97Cdkc1KGRoa2lUwREREVJrvhF1pko8bN27sOiAiIiIqVrFK/MUXXyxalsvlYJomDhw44M7SBQAPHjzA/fv3Sz6HiIiIdqdiws7lcmXvV5fadmFhoSFBERERkVfFhB2LxXwNhnLhwoVdB0RERETFKibsasl6c3MTpmmiv78fhw8frjm5q6oKURRhWRYAQFGUitvbto35+Xn09fUBAMLhMEKhUE3HIiIi2gsqNjq7ffs2nn/+eTz//PO4ceOG29gsm826yxcXF/Grv/qrePXVV2s64OzsLERRhCzLUBQFGxsbSKVSZbe3bRuTk5OYmppyE/v8/Hytr4+IiGhPqFjCjkQisG3bM3Y4AExOTgIA/vzP/9xdZpombty4UbXR2e3btzE1NeU+Hh4eRiKRgCzLJbe/dOmSpwQeiUQgSVLFYxAREe01VUvYhck6m81ieXkZc3NznuWiKCKXy1U8mGEYRcsEQYCu62Wfs7S0BEmSYJomDMOAIAglR1kjIiLayyom7EAgULQsmUxCEAT09vYWrTtw4EDFg1mWhWAw6FlW+Difk+AzmYy77OLFi7Btu+JxiIiI9pqKVeKlxgZPpVJlq6+rlbArJVrbtotGUTNN0/3bKVUPDw/j0qVLuHbtWsVjFcpP+tR5VlZW2h0CleF8D9T7Hn36yScAgP/5F3+BLz/2WMPiIi9eQ3tf1X7Y6+vreOaZZwDsVE/ruo7vf//7RdsuLy9XbbktCILbMtxR+Lhwe2CnVbhDFEUsLS1VPE4p4XAYhw4d8v08ar6VlRUcO3as3WFQGYFAANlstu73aOPdvwcAfOPrXy9Za0e7x2uosz18+LAhhcaKCXt8fBwzMzO4f/++O5LZ3Nycm8CBnUSdTCbx9ttv47XXXqt4sGAwWFTKdh6XGqPcKVXnr3P+LlUiJyIi2quqztYVj8cBoOSUmtlsFr29vYhGo4hGo1UPFgqFipKsZVllW32LoghBEDzJ2fmbyZqIiPaTmif/KNUyOxAIQBRFz79qIpGIp991Op32dNsyTdOzfmJiAslk0n18584dTExM1Bo2ERHRnuB7tq7disfjME0Tuq5D0zT09fV5GrGlUilomuY+jkajsG0bqqpCVVX09PTUVJonIiLaS6pWiTdDpYRbqnqdCZqIiPa7lpewiYiIyD8mbCIioi7AhE1ERNQF6krYS0tLeOGFF3Dq1Cl32RtvvNGwoIiIiMjLd6Oz27dvY2NjA7FYzDN06NmzZ3H37l1PEiciIqLG8J2wg8EgYrEYgOJhRauNJU5ERET18V0lXmmEMc6iRURE1By+E7ZhGFhfXy9avr6+XnEiDyIiIqqf7yrx8fFxTE5OYn193R2a1DRN9PT04MaNG82IkYiIaN+ra6Szubk5mKaJ5eVl2LYNRVEwODjY6NiIiIjoc74T9vLyMgYHB2ue7IOIiIh2z/c97FdffRXb29vNiIWIiIjK8F3C7u3tha7ryOVyEASBVeFEREQt4Dthv/766+7f2WwWS0tLOHDgACRJwuHDhxsaHBEREe3Y1VjigUAAQ0NDeOaZZ/DzP//zePXVVxsVFxEREeXxXcLe3NxEb28vgJ1hShcXF7G9vY3x8XFEIpGGB0hERER1JOzJyUmEw2Ekk0lEIhFcuXIFR48ebUZsRERE9DnfCdu2bYTDYbz22mvNiIeIiIhK8H0POxqN4uzZs82IhYiIiMo4kGvgFFtvvPFGxyXzhw8fIpPJ4Hd/93c9Y53/zM/8DM6fP49/+Id/wLe+9a2i5509exaKouD999/HxMRE0fpvfetb+Nmf/Vncv38fk5OTResnJiZw6tQp/NVf/RV++Zd/uWj9xYsX8VM/9VPIZDL49re/XbT+l37plzAwMIB79+7hN37jN4rWf/vb30Y4HMaf/dmf4dq1a0Xrf/3Xfx1f/epXcffuXSwsLBStn5ubw5NPPok/+qM/wne/+92i9QsLCzhy5Ag0TSs51/l3v/tdfPnLX8atW7fwx3/8x0Xr//AP/xAA8J3vfAd/8id/4ln32GOP4fd///cBAL/1W7+FVCqFQCDgrv/KV74CVVUBAL/2a7+GlZUVz/N/7Md+DL/zO78DAJiZmcHa2ppn/Y//+I/j6tWrAIBXXnkFf/3Xf+1Zf/ToUcTjcQDAyy+/jL/927/1rD927Bh+5Vd+BcDOD9QHDx541h8/fhy/+Iu/CAD4uZ/7OXz44Yee9SdPnsQv/MIvAADOnDlTdG667bNnGAY++eQTfOMb3wDg/7P3f/7u/+HRL/fg8OM5fOlLX+qoz146nfas79bPXjabRSAQ2HOfvULd+r0XDAbx0ksvIRwO49ChQ0XPq1XVKvFXX30Vo6OjeOaZZwAAL774YsntcrkcDMPouIRNRES0F1QtYScSCfz0T/+0m7BfeOEFxGIxT4kof9u5ubnmRFonp4S921821DwrKys4duxYu8OgMs6cOeOOuVCP7/9BEj/61X+Dn/yJXMnvDdo9XkOdrVF5qGoJOxaLeR6/9tprZccQL9yWiIiIGsN3o7NKE35wMhAiIqLm8N2ta3t72zMEqWmaWFtb47jiRERETeS7hK1pmuexKIoYGhrC4OAg7t6927DAiIiI6Au7Gku8UH63KSIiImqcqlXi2WwWyWQS6XQa29vbME0Tuq4XbWeaJhRFaUqQRERE+13VhB0IBHDu3DmcO3cOs7OzePLJJzE6Olq0nSiK7LJBRETUJL4anU1NTWFpaaloso9sNgvLspiwiYiImsT3PeyhoaGiZYFAAMFgkI3OiIiImsR3ty5gZ07s5eVl2LbtWb6xsYFTp041JDAiIiL6gu+Evby8DFVVcfToUdi2DUEQAOxMu/nKK680PEAiIiKqI2Gvra3h9ddfB7DTMhz4YoSz9fV1d8xxIiIiahzf97B7e3vdv0VRxNtvv+0+3traakhQRERE5FX3wClOA7PV1VVsb28DQNHcsERERNQYdbUSv379Ou7cuQNgZ5L15557Dt/85jdZwiYiImqSulqJj4+Pu3+Looh33nkHa2trRf2ziYiIqDEaNpb40aNHcePGjUbtjoiIiPJULGG/+OKLNe8ol8vBMAxfzyEiIqLaVEzYuVwOsVis5iFHE4lEQ4IiIiIir4oJOxaL+bovHYvFdh0QERERFat4D9tvI7L19fWatlNVFalUCpqmQdM0X8eYmZnxtT0REdFe4LuV+BtvvFFyuW3b0DSt6ljis7Oz6O/vhyzL7uNUKuU+rvZcZ3Q1IiKi/cR3wl5YWIAkSZ772rZtY3NzE4qiVH3+7du3MTU15T4eHh5GIpGomrANw/AbKhE10fvvv4/79+/jo48+qrhd39P/GMBf4t13WxPXfrWystLuEPalgwcP4sknn8SRI0eafizfCTsajeLcuXMl1y0tLVV8bqmkKwgCdF2vetxMJoPjx49zNDWiDvD+++/DNE089dRTePzxx/HIIw3rIUrUNT777DN88MEHeO+99wCg6Unb91VWLlkDcGfuKseyLASDQc+ywselpFIpRCKR2gIkoqa7f/8+nnrqKRw+fJjJmvatRx55BIcPH8ZTTz2F+/fvN/14dY10Vs7m5mbF9YXzZxeuK5XwneXVfgxUk8lkdvV8ai5W53WubDYLoPg9evzxx9sRDlHHefzxx/HRRx81/XvMd8IuNTDKgwcPsL29XfUetiAIsCzLs6zwcaFkMlnTvfFqwuEwDh06tOv9UOOtrKzg2LFj7Q6DyggEAshms573aGVlhSVros8510K577GHDx82pNDoO2E/ePAAV65c8SwLBALo6empOsBKMBgsKmU7j0uVoA3DgCRJfkMkIiLac3wn7CtXrpTtn725uemZL7tQKBQqSsyWZZVNypZleRqkra6uwjRNqKoKWZYhiqLf8ImIiLqS74TtJGtnDux8qqritddeq/j8SCTi6XedTqc9Vd6macIwDMiyDEmSPMlc0zRsbm4iGo36DZuIiKir+U7YS0tLmJ6eRk9PD3K5HADgwIEDyOVy2NzcrJqw4/E4VFWFruswTRN9fX2ePtipVAq6rhf1y9Y0DalUyi1hK4qy64ZoRESllGsES9ROvhO2aZp45513Sq67fv16TfuoVEKORqMl1yuK0pDGZ0RE1Vy6dAnDw8M1jcCYzxk6OR6PNyMs2ud8J+xQKFR23fj4+K6CISLqBMvLy8hms74Ttt/tifxoaL+M5eXlRu6OiKjlUqkUYrEYdF2vOHZEKYXtbogayXcJe3BwEL/5m78JAOjp6XHv89i2jTt37uB73/teYyMkImoh0zQRjUaRSCQaNg4EUSP4TtiJRAKmaUIURTx48AAPHjxw1zkjIhHR/nQ79af4gztvtTWG/zB8Aufk5+p6bn5jM6dHS6mEraoqQqEQbNuGruuIRqOwbRuJRAIAcPPmTQA7Y0lYlgXbtrG6uorjx4+7JXDDMJBIJBAIBDA6Ogpgp9dMX18ffyRQSb4Tdn9/P2KxWNl1RETdKplMuvMWyLKMsbGxohbjmqZBFEVP1bdt2wiFQohGo1BV1V0+PT2NiYkJyLIMWZYxMDCAe/fuAdhpD6QoChKJBERRhCiKCIfDGBgYYMKmknwn7EpdHYaGhnYVDBF1t3Pyc3WXbjtBfnKWJAmCIBRVi4uiiJmZGXfQp0oNzebm5ooGeMo/hvO/s03+LUZ2K6NCvhud9fT0lJ3k48aNG7sOiIioHUzTdMd5cP6Fw2FomubZTpIkxONxpNNpjIyMYGRkpGzjtGAwCFVVoWmaO2pj4fwJHLGRauW7hP2d73wHm5ubyGaz7hjiwM4Y4/fv3y85OQgRUafTdb2o/7QkSRgZGXHb7Tjb5bcGn5mZgaZpJcePGBkZwdzcXFF3WJagqR6+E3Y2m0UsFiuaxzqXy2FhYaFhgRERtdLGxkbRslAoBFEUkUql3IRsGAYAuAlbURTPnAcOp8GZk6zzS+H5ozlWm7GQyOE7YcdisbKTf1y4cGHXARERtZKu61BVFZlMBv39/Z570pqmwbIstzDiDIlsmiZSqRSAL7qBGYaBxcVFZDIZaJoGRVEQiUTcFuUAcPnyZXdoZcMwoKoqTNOEpmmIRCKYn58HsNMbJxqNsrqcPA7knAHBG+Du3bs4depUo3bXEM48pJwPu3NxPuzOdubMGWSzWSwtLbnL+J4ReVW6JhqVh3yXsN94442Sy23bhqZpHZewiYiI9gLfCXthYQGSJCEQCLjLbNvG5uYm+w4SERE1ie+EHY1Gce7cuZLr8qvMiIiIqHF898Mul6yByoOqEBERUf0aOltXuQFViIiIaHd8V4mXGhjlwYMH2N7e5j1sIiKiJvGdsB88eIArV654ljkjnuU3RCMiIqLG8Z2wr1y5UnbgFCIiImqOigl7e3u7aBmTNRERUetVbHSWTqdx+vRpDAwMYHZ2tmi8XNM0sbS0hLt37zY1SCIiov2uYgl7aGgIlmUhEomUvD/tTLoO7IyAdvbs2eZESUREtM9VLGHfvXsXw8PDNTUmO3v2bNlhS4mIiGh3KiZsy7Jw+PDhmnfWwHlEiIiIKE/FhJ0/f2ststnsroIhImq1VCqFkZERfO1rX8Ps7Kz7vadpGk6ePOm24XGYpomxsTGcPHkSmqbBMAyMjY1hbGys4nFmZmYwMzPT0NhN08Ts7CwGBgaKjm8YBi5evFgUfzOYpglVVTEwMICRkRGoqlpxOdCc87EbnRZPSbkKEolEpdVFrl+/7mv7Vvjwww9zP/jBD3Iffvhhu0OhMn7wgx+0OwSq4Pnnn8+dOnXKs2yvvWeZTCb39NNP5yzL8ixPJpO5Z599tmj7jY2NXDqd9mx3/vz5isdIp9Oe5+RyudzCwsIuovbu59lnn80tLi6WXNcqJ06cKHm806dP565evepZVup8tNNu46l0TTQqD1UsYQeDwZpbgC8vL7NKnIi6UigUgiAISCaTnuWCIMC27aIeMoZhQJIkz3bVSJLkeQ4ArK6u7iJqb5yXL1/GzMxMUc2o0zC43Xp6ejyPS52Pduq0eEqpmLDHx8exuLiIt99+u+JO1tfXcf36dYyPjzc0OCKiVolEIkilUp5ltm1DkqSi5Y3Q6GpqWZYhSRImJycbul/qHFVHOnvttdcwOTmJr3zlK5BlGeFw2F2XyWSQSqWwtbWF3/7t325mnETUBf7rf/vP+NOV/9jWGJ479u/xb//lv/P9PFmWMTY2Btu2PSVmWZaRSCQQj8cB7CTxcqVWXddh2zZWV1fR19fnzq9gGAYSiQQA4ObNm9B1HaZpYm1tDaqqQhAEd1vbtjE/P4/+/n6srq7i+PHjNZf85ubmcOLECaRSKciyXHIb27ahaZr7GkzTRDQa9cQZCAQwOjoKYGc8jvzX0giF5wMAVFVFKBRyazSi0ShM00QikUBvby+OHz+OYDBYdG6d/VmW5Z77/HNWy2sqFU/+++DUtJQ7p61SNWGLoog333wTqqpiYWEBpml61imKwpI1EXU9SZLcanFFUaDruvulPzMz4z7Wdb3kF3cmk3HHppAkCQMDA25CCIVCiEajbqMrSZLcqmsnWTpGRkbw5ptvQhAEyLKMkZER3Lp1q6Zqd0EQEIvFMD097b6eQufPn/fsL5VKYWZmBvF4HKFQCIqiIJFIuK8lHA57Xks1hbcPAHjyRqnz4fyAyP9h4tRuTExMIJFI4Nq1awB2fkBdvHgRANyYpqenMTExAVmWIcsyBgYGcO/ePfdY1V5TYTzOebp8+TJCoRAA4OTJkxBF0X3cDjWPJR6NRt0PlmmaHXNfhIg6x7/9l/+urtJtp3CqxRVF8ZS0nWrxSiXdYDDofi86zyssrVeTSqUgCILnOeFw2P0RUQtFUZBKpXDp0iU3yeXvPz8+YCcBTk5OIhaLeY5d72uRJKnoR0hh24BCoihiZmYGlmVBkqSiH0SFQ2KPjo5icnLSPSdzc3NFOSk/Xr+vyakByU/Oc3NzbU3WQJ3zYTNZE9FeJMuyW61duDyZTFZMWo34XnRKorquu/+ce9N+xONxLC0tFd17N00TwWCwaHtBEIpqT8txurQ5/0qVqP2SJAnxeBzpdBojIyMYGRmp2K1YFEXP+mAwCFVVoWmaG49lWUXPqVWpQmm7kzVQx2xdRER7lVONfOnSJc80wpFIBDMzM5572Y3kJEsnEe22tbIoim7VeCwW8ywvTGRA5fvyhZx7vI3k3G7IvwWhaVpRSd1RmFBHRkZKloD91nA4RFEsqsbvBHWVsImI9qrBwUFks1nPF70gCFWTaKlEWIkoiu5gU04CkmW5qLRrmmbVUuzGxkbRsmg0WlQSlWUZtm179p9KpTA0NOR5vX5fSy22trbKrjMMw/MaC6v/19bWPI9VVfU0lLMsy03W+a83f59+XpMkSRBF0fN8wzBgGEbN+2gGlrCJiPKMjo6WLF3JslyyFGoYBlRVhWma0DQNkUgE8/PzAIBEIoFoNArbtrG4uIhMJgNN06AoCkKhEERRhKZpnv3dunXLbZ0MVP6x4Ix0try8DACYmpryrJ+bmyuqFn/zzTfd/Tstq5173bW8llLnwDRNpFIpmKaJZDLptnrPXw58kWgNw/CcD+dHihNrfst1AO55CgaDME0TiqK497lDoRAikYjbyhwALl++DFVVoShK3e/PrVu3kEgk3NsItfxoa7YDuT0+2snDhw+RyWQQDodx6NChdodDJaysrODYsWPtDoPKOHPmDLLZLJaWltxlfM+oVVKpFO7cuVPUgK7TVLomGpWHWCVOREQdjfNU7GDCJiKijqTrOjRNc6uq9zvewyYioo7UDeN7txJL2ERERF2ACZuIiKgLMGETERF1ASZsIiKiLtCWRmeqqnqGyKs0qL0zFRwAd9q0Rk7zRkRE1A1anrBnZ2fR39/vjlIzOztbce7W+fl5z+g9J0+eBFA5yRMREe01La8Sv337tic5Dw8Pl+1fVzjmLbCTqPPnLCUiItoPWpqwSw2cLghCxYHtl5eXPUm7cGB8IiKi/aClVeKWZRXNxVpqblaHIAi4d++eZ1k6nWZHeiIi2ndamrArTUhey7yltm1jeXkZt27d8n3sTCbj+znUOisrK+0OgcpwxnHme0RUWbOvkZYmbEEQiuYk9TNH6aVLl0pOUl4LztbVuTjzU2cLBALIZrOe92gvJe9UKoWFhQUYhoHx8XFcuHABgiBA0zSoqgrLsnDu3Dm38atpmpiZmXGngAyHw0gkEgCAmzdvlj3OzMwMACAejzcsdmdKyoWFBYiiiEgkgmg0WnZ5s+LYjU6LZzeqzda1Wy1N2MFgsKiU7TyuVrpWVRWjo6OsDieihnLmuR4ZGXGTNbDTwDUYDGJ6etrTU0UURcTjcZim6X4fKYpSdXKKUj1hnPmh6yWKIqLRqDvPs7MvZ3kymcTg4KDnGOV65LRLp8XTyVra6CwUChUlZsuyqibhVCqFUCjkblepkRoRkV/Od1MymfQsFwQBtm0XfecYhuH53qpW4ABKT2Sxurq6i6hr09PTUzWOduq0eDpZy7t1RSIRpFIp93E6nfb0qXaqchy6rsOyLITDYbebV6nW5kREu1H43QTs1ABKklS0vBFmZ2cbvk/a21o+cEo8HoeqqtB1HaZpoq+vz1MlkkqloOs6ZFmGbdsYGxsD8MV9DgAYGhpqddhEVIP3/8f/wvv/fa2tMRz5yaM48i9+wvfzZFnG2NhYUQNYWZaRSCTce6y2bUMUxZL70HUdtm1jdXUVfX19bmHEMAzPfW7n+29tbQ2qqkIQBHdb27YxPz+P/v5+d3THRpVAC+MAdqrlQ6GQW5Pg3ANPJBLo7e3F8ePHEQwGi16Tsz/LstzXnB+rc6xAIIDR0VEAOwW0Suel8PU7NRysNt/RlqFJK92ziUaj7npBEPCXf/mXrQqLiPYxSZLcanFFUaDrupt8ZmZm3MdOgaJQJpOBKIoQRRGSJGFgYMBNTKFQCNFo1B30SZIkt/1O4ffhyMgI3nzzTQiCAFmWMTIyglu3blWtdi91q7BwzIrCODRNc+N1OLUKExMTSCQSuHbtGoCdHy4XL14E8MVIk9PT05iYmIAsy5BlGQMDA25X3FAoBEVRkEgk3PMSDocrnhcAOH/+PC5fvuw2Lj558iREUayrsfFe05aETUR705F/8RN1lW47hVMtriiKp6TtVItXKukGg0G35O08r5buqvlSqRQEQfA8JxwOuz8iKpEkqSj5F96TLySKImZmZty2RIU/RI4ePep5PDo6isnJSTeWubm5otqG/Nfs/F/reXFqHvKTc709g/YiztZFRPQ5WZbdau3C5clksmICLldN7odTItZ13f0ny3LTGmVJkoR4PI50Oo2RkRGMjIxUHC9DFEXP+mAwCFVVoWmaW8Iv7Krr57yYplm0PZP1F5iwiYg+51SLX7p0yZMkI5EIbNtGIpFoSvI0TdNNVk6VdP6/RvwYKMWp5r927Rru3buHcDhcsXtaYUIdGRmBJElQFKWoWr0eoihy6OkKmLCJiPIMDg4im816StKCIFRN1H4GgQJ2kpMzipyTCGVZLpovwTTNXXVl3draKrvOMAzPvgur3dfWvA0I8/uNOw3OnBJwfpLO36ef8+L8OMl/vmEY7Bn0Od7DJiLKMzo6WrKU5wywUsgwDKiqCtM03QFM5ufnAQCJRALRaBS2bWNxcRGZTAaapkFRFIRCIYiiWFSivXXrlttKGqj8Y8HpBmuaJpLJpNvaPH858EWiNQzDE4fz48DptuaM3uZw4gsGgzBNE4qiuPe5Q6EQIpGI28ocAC5fvgxVVaEoSt3n5datW0gkEjBNE8FgsKYfS/vFgVwul2t3EM3kDAnHoUk7F4cm7WxnzpxBNpvF0tKSu4zv2d6XSqVw584dt5U4VVbpmmhUHmKVOBERleRU2VNnYMImIiIPXdehaZpbVU2dgfewiYjIg+N7dyaWsImIiLoAEzYREVEXYMImIiLqAkzYREREXYAJm4iIqAswYRMREXUBJmwiIqIuwIRNRETUBZiwiYiIugATNhERURfg0KREtK+ZponFxUVcv37dnTIS2JlH2jRNjI6OusN0GoaBRCIBALh58+aujqvrOhKJBHp7ezkjFtWECZuI9jVRFDE1NYWlpSUoigJFUdx1tm1jYGAAN2/ehCRJCIVCiEajUFV118eVJAmKoiCdTu96X7Q/sEqciAiAIAgll0mShMXFxaYcMxgMNmW/tDcxYRMRVZDJZNDf39/uMIhYJU5EjfN3//dj/O//+3FbY/in/+hR/JN/9Oiu92OaJlRVxcTEBKLRaNntDMOAZVmwbRurq6s4fvy4Z2pK27YxPz+P/v5+CIIA27Yhy3LRfmzbxokTJxAOhxGNRjm9JRVhwiYi+pxhGEilUu5jURRLVpXnm56exsTEBGRZhizLGBgYwL1799z158+fx+XLlxEKhQAAJ0+ehCiK7mOHZVmIxWKee+hE+Ziwiahh/kmDSrftEgqFikq/Y2NjMAwD8Xi85HPm5uYgiqJnmW3bEAQBuq7DNE1Pcp6bmytK1oZhQNf1iiV5It7DJiKqQFEUaJpWdn0wGISqqtA0DbquA9gpLQM71eqFybwwWa+trSGTyWBhYQG2bTc4etpLWMImIqrASb7ljIyMlCw127YNURRhmmbF5zvdu2zbxqVLl9gnm8piCZuICChbutU0rex9ZafBmZOs8/eh6zokSYIoim7J23mOYRhF+4pGo1hbW/PcQyfKxxI2Ee1rpmkilUq5/ztJd2trC2traxgcHMTU1BSAnWS7uLiITCbjJvJIJAJVVd2kffnyZaiq6ib5W7duIZFIwDRNBINBt2+3ruvQNM09rizLEAQBk5OTiMVivJ9NRQ7kcrlcu4NopocPHyKTySAcDuPQoUPtDodKWFlZwbFjx9odBpVx5swZZLNZLC0tucv4nhF5VbomGpWHWCVORETUBZiwiYiIugATNhERURdgwiYiIuoCTNhERERdgAmbiOry2WeftTsEoo7QqmuBCZuIfDt48CA++OCDdodB1BE++OADHDx4sOnHYcImIt+efPJJvPfee9je3mZJm/atzz77DNvb23jvvffw5JNPNv14HOmMiHw7cuQIAOBv/uZv8NFHH7U5GqL2OXjwIERRdK+JZmLCJqK6HDlypKYvqe//QRI/+tV/g5/8iRwCgUALItt/OPLc/sAqcSIioi7QlhK2qqoQRdGdtq7cTDj1bk9ERLTXtLyEPTs7C1EUIcsyFEXBxsZGxenk/G5PRES0F7U8Yd++fRuyLLuPh4eHoWlaw7Ynova7873/jDt/8KdY0v4MR3706XaHQ7QntLRKvNSk7YIgeCZ33832RNQZvvSlf4bHxKeR+/hD5D79BB/+/Xv48k+G2h0WUVdracK2LAvBYNCzrPDxbrYvxZnum11POtvDhw/bHQKVEQwG8SM/8iO+3qNHvgTktv8O//q4+PmSI/j000/x6aefNidI4jXUwZz84+SjerU0Ydu2XXGdIAi72r6Ujz/+GADw7rvv1hgltUMmk2l3CFTGSy+9BMDfe/SVfwoAFjIZqzlBURFeQ53v448/xmOPPVb381uasAVBcFt6Owof72b7Up544gk8/fTTePTRR3HgwAFfzyUiItqtXC6Hjz/+GE888cSu9tPShB0MBotKzc7jUqVlv9uX8sgjj3CwBiIiaqvdlKwdLW0lHgqFihKtZVmQJKkh2xMREe1VLe/WFYlEPP2o0+m0ZyAU0zQ966ttT0REtB8cyO222VodVFVFKBSCaZoAvCOXqaoKXddx8+bNmrYnIiLaD9qSsImIiMgfTv5BRETUBZiwiYiIugATNhERURdoy/SazWbbNjRNw9bWFqampqpuz+k7W8/POdc0DYZhuJPApFIpRKNRiKJY9jnkD6e87Wy8Xjpbq3LOnith67oOXddhmiay2WzV7Tl9Z+vVc86TySTGxsaQSCSgKAq/fBqIU952Nl4vna2lOSe3R129ejU3PT1ddbtnn33W8ziTyeTOnz/frLAo5/+cLy4uNjukfc3v+8FrprV4vXSHVuScPVfC9oPTd7Yez3ln4ZS3nY3ne2/Z7fu5J+9h16oR03eSP/Wec03TEAwGec+0wdox5S3VjtfL3rLb62dfJ+xGTN9J/tRzzsPhMARBcO/DXbx4EcFg0G1UQ/Vrx5S3VDteL3vLbq+ffV0l3ojpO8mfes55KBTyNJrp7+/HwsJCU+Lbb9ox5S3VjtfL3rLb66fjS9ipVAp37typuE1PTw/i8bjvfTdi+k7y9x7Vc851XffM0CaKYsl7QeRfO6a8pdrxetlbdnv9dHzClmW5aVU5nL6zMfy8R37PuWmaGBsbw7179zzPYzeVxuCUt52N18vestvrZ99ViXP6zvbzM8WqKIqIxWKeD/mdO3cQjUZbF/AexylvOxuvl+7WyOtnz83WZRgGdF2HpmkAdlpHSpKEUCgEgNN3dgo/U6zmf+C3trbQ19fH96jBOOVtZ+P10rlamXP2XMImIiLai/ZdlTgREVE3YsImIiLqAkzYREREXYAJm4iIqAswYRMREXWBjh84pZOpqoqFhQUEg0G3Wf7W1hay2SxCoVBR14ty2wYCgZomPXfkd9vI7z7QbqVeI7ATr6IoLYlT0zSoqopYLNa2sZNVVYWmaRAEAZFIxF1umiZ0XYcgCHjzzTfbEluz8Frwyv8MFL7XzmfU6TPdiJgNw0AikQAAT/ehbuCcK2Cne5MsyxzopRzfk36Sx8svv5y7evVq0fLp6emiuVHLbZtMJn3NJ3z69OmcZVm5hYWF3MLCgv+gm6jcazx9+nQumUy2JIbp6emWHauccufBsqy2zR1tWVZT989rwWtxcTF34sSJkvNTLy4uNvz9KJxXudb5mTtBuc9DI5U6H82+JhqNVeIN0NPTU7QsFotB0zS3Y3ylbZ2hPWdmZqoey7Ztd1aXaDTakSMYlXqNExMTmJ6ebsnxA4FAS45TD0EQ2jaMZzKZbPoxeC18IRgMIh6PuyXffKIoNnzs9cL9DQ8Pc8CUPKXORyuuiUZiwm4S5+KpddB9RVGQTCaLvtQKWZbVldVFgiC4X7D7XSgUast5SKfTLT8msL+vBUmSEA6HMTs72/Jjh0KhjrlF0AlKnY92XRP1YsJuEucL2c8FEw6HPWPM7iXpdBqSJO3rGZ2chBUOh1t+7NnZWWSz2ZYfF+C1EI/Hcf369ao/QKi12nlN1IuNzppkfn4eQ0NDvkoAoihidXW17Hpd15FOp2GapttoRRRFTE9PIxwOQ5ZlmKaJdDqNa9euAdhp4BIMBgHslEicKiGnkUogEMCFCxdgWRZs28bq6iqmpqag67q7nSiKdTfgsm0bmqYhm82WbAzjfCmbpolQKORWF1eLL5/zGp1SfKmLsN3nwTRNd97b/B8tleKq9L46YxHbtu15nqZp7mfOWRcMBpHNZt3PDbBTim3Vj6f9fi2IoghFUTAzM1OxQVi98ZVimqZ7W8E5pp/9lPt8AZWv2UrnvxHKnaPC9U5sTrzRaNRzPlKpVFuvibq1+yZ6t3v55ZdzL7/8ci6ZTOaSyWRucXExd/Xq1ZINTV5++eWKDWOuXr2aO336dMXjbWxsFDXKSafTudOnT+c2NjbcBjjO8TY2NtztChs85T8vP8bCxh/PPvtsxZjyOecjnU7nkslkbmFhIXf16tWSjTsKz8X58+c9sZSLL51Ou49LnevCBm7tOg/nz593X/+JEyc8cfuNq/B9zWQynnPgfPYK9+e8hsIGSc3Aa8Er/zNoWVbu2WefdT8Du/ks5D8nfz+F56PUe17Lfsp9vnI5f9ds/vmvppZGZ9XO0fT0tPs6LMvKnThxwvP8wvPRimui0Vgl3gDOr25ZlqEoCqampupu7NHb2+v7Oc6k6E5Dlmg0CsMwsLm56SnVCIIAURTdLhT5z8t/LeX2X6v+/n5IkgRZlhGNRjE8PIyRkZGi7VZXVz3VnkePHnVLM5Xic6oWbdvG7du3i851fpVzO8/D0aNHEY1GS34e/MblvK+maWJtbc1TvTw8POw+Z3Fx0bO/4eHhmuNtBF4LpQmCgFgsVrIx3W7iq1TNXqq0WG0/1T5ffq5Z5/w3Qi3nKJlMutd+qXYTHV96rgGrxDuIaZp1N6IpfF4mkyn5hSeKYtUPcanWu7sRCoVgWRZ0Xfe0kM6vKjNNE9lstujLsNJFput61fPVKedBURRkMhnfcRW+Pqcfd/6XpG3bbl/3ixcv4mtf+5r7g6lbWwnvxWtBURS3z3F+jLuJrx7Vrqlyny+gtmu20Q0BdV2HaZpVz5EoirAsy3193dAo0S8m7A6yvLyMCxcu1PXcwq5MlUoBW1tb7t/O/aB8zfgl6txXyufcP3Lug5XqjlUqPj865TwUdueqNa5S72tvb29R1zBZlmHbNq5duwbbtpHJZKCqKgzDQDweLzpGtYQ4Njbm3nMHgEgk0tJuU3v1WojH45icnEQsFmtIfPWotJ9Kny+gtmu20d0qq/Uucc6RoihYXFzEhQsXkEwmMTEx4ev9282PxFZhwu4QqVQKg4ODDeuGIUlSyT6Gpmni+PHjDTmGH4IgYGNjw7NsZGQEb731VsmLqpZqx/wJ4MvptPPgqDeuUChUtu/o/Pw8pqam3B8HkiRhbGys5LZOA6py2jla1l6+FpxuXqqqurUfnRRfpc8XsPtr1i+ndF/LOQoGg7hw4QIymQwikYjvH1vVrolOwHvYLZb/i9lhGAY0TcOVK1fq3m9hy+hQKITe3l5PlZpT8sqvJs0vReVvtxulXmM4HMba2hqAndf71ltvAfCWYJzX4LRMrRafKIqIRCKe6jtgpwrNeV47z0MltcZV+L5KkoRgMFjUp9m5j+f873C+gArvd3bC/by9fi2Ua+Uej8c970Wr46u0n0qfL2dZtWu2kV2lLl26BEEQajpHzvdArd1HO/GaqIYl7F1QVRVra2tYW1tDT09PxSpDVVWxvLyMzc1NtxtBvlpKNKZpYnZ21q3uVBQFpmlifn7es8z54F27dg2qqnoalNy6dQvAzhfj/Py8W8UVjUaRSqWwtLTkXiCSJLnPTyQSiEajFX+BOoND5O/TEYvFcOnSJbfBiizLOHfunFu9FgwGEYvFkEgk3HvTtcQXj8ehqqo74pVt25Akyb1PKElSy89D4eei3NjItcRV6n29efMmVFVFJpNxqzed+6OiKLpfXLZtu++BIAhQFMXt9tLocdZ5LXjNzMy43RljsZgnGYiiiPHxcc/2u41PFEXP+ZBl2X2saRoURan5dZb7fAGo6Zotdf4rcT4PzmfHef3OvWun4WSlcwTsfKeMjIy43wMAMDQ0hKmpKc/nxTkfzb4mmuFALpfLtTsIIiKiehmGAV3Xi34oqqqKra0tXxPKdDJWiRMRUVfTNK1kCTkajbq34vYCJmwiIupqkiQVtWUBUNSVtNuxSpyIiLqec887f2heAF1xb7pWTNhERERdgFXiREREXYAJm4iIqAswYRMREXUBJmwiIqIuwIRNRETUBZiwiYiIusD/B68qN9ACy+C1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 540x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfQAAAFRCAYAAACCB1/XAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA4M0lEQVR4nO3df2wbd343+Ld313FahzPadHHPZaMxDote9olJtc89jrDwGOhTxIY9VFu0VmyPD+gWViLKRQBbKEClLWopG9pAuxaLVi5QrDR2bHSL5zR2N73iKUzK2PR5ngIcZePVcweII1+D7u4DjdK99h7EmqGSjde74f2hzoRDUiSH4i+N3i8giDkzHH7mqxl++P3Od77fPcVisQgiIiLa0T7T7QCIiIho+5jQiYiIQoAJnYiIKASY0ImIiEKACZ2IiCgEmNCJiIhC4HOd/kDHcaDrOtbX1zExMVF3e03TIEkSbNsGAKiq2u4QiYiIdpyO1tANw4BhGLAsC4VCoe7209PTkCQJiqJAVVWsrq4im812IFIiIqKdpaMJXZZlKIqCSCTS0Pa3b9+Goije66GhIei63q7wiIiIdqyevYdummbFMkEQYBhGF6IhIiLqbR2/h94o27YhiqJvWfnrRnzyySf48MMPsXfvXuzZs6dV4RERETWkWCzi8ePH2L9/Pz7zmfbVo3s2oTuOU3OdIAgN7efDDz/Ee++916qwiIiImvLcc881fMu5GT2b0AVB8Hq2u8pfN2Lv3r0ANgvyiSeeaElsYZDP5xGLxbodRs9huVR69dVXsbGxgb/4i78I/F4j9w4+F/kFrK8t4vjQ0TZE1108XyqxTCr9+Mc/xnvvveflo3bp2YQuimJFLd193WjtHIDXzP7EE09g3759rQswBFge1bFc/GzbRqFQaKpcip8UUcRe/PQnPw1tuYb1uLaDZVJdu2/79mynuGg0WpG4bduGLMtdioiIiKh39VRCtyzL95x5PB73vc7lchxYhoiIqIqOJnTTNKFpGhYWFmAYBjRN8z2els1mfc+Zp1IpWJYFwzCg6zoOHDjgey6diIiINnX0Hno0GkU0GkUikai6PpFIVKzbalsiIiL6VE81uRMREVFzmNCJiIhCgAmdiIgoBJjQiYiIQoAJnYiIKASY0ImIiEKACZ2IiCgEmNCJiIhCgAmdiIgoBJjQiYiIQoAJnYiIKASY0ImIiEKACZ2IiCgEmNCJiIhCgAmdiIgoBJjQiYiIQoAJnYiIKASY0ImIiEKACZ2IiCgEmNCJiIhCgAmdiIgoBJjQiYiIQoAJnYiIKASY0ImIiEKACZ2IiCgEmNCJiIhCgAmdiIgoBJjQiYiIQoAJnYiIKASY0ImIiEKACZ2IiCgEmNCJiIhCgAmdiIgoBJjQiYiIQoAJnYiIKASY0ImIiEKACZ2IiCgEmNCJiIhCgAmdiIgoBJjQiYiIQoAJnYiIKASY0ImIiEKACZ2IiCgEPteND9U0DZIkwbZtAICqqnW3FwQBAOA4DhKJRNtjJCIi2kk6XkOfnp6GJElQFAWqqmJ1dRXZbHbL7TVNQyKRgKqqUFUVsixD07QORkxERNT7Op7Qb9++DUVRvNdDQ0PQdX3L7TOZjO91NBrF8vJy2+IjIiLaiTqa0E3TrFgmCAIMw9jyPaIo4uLFi95rXdcxNDTUlviIiIh2qo4mdNu2IYqib1n563KpVAqLi4sYHByEpmkQRdFXwyciIqIOd4pzHKfmOrfjWylJkjA2NgbDMJBOpzE6OtpUQs/n84HfE3ZLS0vdDqEnsVz8CoUCgObK5V/+5Z/xRRH48MMPQ1uuYT2u7WCZdEdHE7ogCF7Pdlf563JTU1NQVRWJRAKGYWB8fByWZeHatWuBPjsWi2Hfvn2BYw6rpaUlHDp0qNth9ByWS6VIJIJCodBUudgPN38M7N+/P5TlyvOlEsuk0qNHjzpSqexok7soihW1dPd1tdq5aZqIRCKIRqMAAFmW8fbbb2NxcbH9wRIREe0gHU3o0Wi0InHbtg1Zlqtub9s2+vr6fMsEQcDhw4fbFSIREdGO1PHH1uLxuO+581wu5xtYxrIsb70syxU94B3HgSRJnQmWiIhoh+j4SHGpVAqapsEwDFiWhQMHDvg6uWWzWRiG4S1LpVKYnp7GgQMHvG0mJiY6HTYREVFP68rQr7WGbk0kEr71kiQxgRMREdXByVmIiIhCgAmdiIgoBJjQiYiIQoAJnYiIKASY0ImIiEKACZ2IiCgEmNCJiIhCgAmdiIgoBJjQiYiIQoAJnYiIKASY0ImIiEKACZ2IiCgEmkroCwsLePnll3H8+HFv2Z07d1oWFBEREQUTeLa127dvY3V1FclkEpZlectPnz6Ne/fu+ZI8ERERdUbghC6KIpLJJADAtm3fumKx2JqoiIiIKJDATe6CIGy5znGcbQVDREREzQmc0E3TxIMHDyqWP3jwoKLGTkRERJ0RuMl9dHQU4+PjePDgASKRCCRJgmVZ6Ovrw40bN9oRIxEREdUROKEDwMzMDCzLwuLiIhzHgaqqOHz4cKtjIyIiogYFTujj4+P4hV/4BbzyyiuQJKkdMREREVFAgRP6kSNHEI/Hq67b2NjAU089te2giIiIKJjAneIkSdqy85uu69sOiIiIiIILXEPPZDIwTROFQgGSJCESiQDYfGTNNE288sorLQ+SiIiIaguc0PP5PCYmJiCKYsW62dnZlgRFREREwQRO6FeuXMHBgwd9ywqFAtbX13H+/PmWBUZERESNC3wPvTyZA0AkEkFfXx/W1tZaEhQREREF09Rz6Gtra94z6KVWV1c5OQsREVEXBE7oi4uL0DQNBw8ehOM43tjujuPgtddea3mAREREVF/ghL6ysoI333wTALzpU90BZh48eIDnn3++heERERFRIwLfQ+/v7/f+LUkS3nnnHe/1+vp6S4IiIiKiYAIndNe9e/cAAMvLy9jY2ACwWXsnIiKizguc0E+cOIHr16/j7t27AIBEIoEXX3wRX/nKV1hDJyIi6pKmermPjo56/5YkCe+++y5WVlaqPtJGRERE7dd0k3u5gwcPYnFxsVW7IyIiogCaqqEvLi5WDCLjOA7u3r2Lb33rWy0JjIiIiBoXOKGn02lYllV1LvRCodCSoIiIiCiYwAl9YGAAyWRyy3VERETUeYHvobsjw1Vz4sSJbQVDREREzQmc0GtNwnLjxo1tB0RERETB1Wxyf+WVVyqWFYtFWJaFPXv2eLOsAcDDhw/x/vvvV30PERERtVfNhF4sFre8X15t27m5uZYERURERMHUTOjJZDLQYDHnz5/fdkBEREQUXM2EXi+Zr62twbIsDAwM4Kmnnmo4+WuaBkmSYNs2AEBV1ZrbO46D2dlZHDhwAAAQi8UQjUYb+iwiIqLdoGanuNu3b+Oll17CSy+9hBs3bnid4QqFgrd8fn4ef/AHf4DXX3+9oQ+cnp6GJElQFAWqqmJ1dRXZbHbL7R3Hwfj4OCYmJrzEPzs72+jxERER7Qo1a+jxeByO4/jGbgeA8fFxAMB3vvMdb5llWbhx40bdTnG3b9/GxMSE93poaAjpdBqKolTd/tKlS74afDwehyzLNT+DiIhot6lbQy9P5oVCAYuLi5iZmfEtlyQJxWKx5oeZplmxTBAEGIax5XsWFhYgyzIsy4JpmhAEoeoodURERLtZzYQeiUQqlmUyGQiCgP7+/op1e/bsqflhtm1DFEXfsvLXpdwfAPl83lt28eJFOI5T83OIiIh2m5pN7tXGZs9ms1s2j9eroddKxI7jVIxCZ1mW92+3Vj40NIRLly7h2rVrNT+rXOmPAtq0tLTU7RB6EsvFz/0eaKZc/uVf/hlfFIEPP/wwtOUa1uPaDpZJd9R9Dv3Bgwd4/vnnAWw2fxuGgb/+67+u2HZxcbFuz3NBELye7a7y1+XbA5u92l2SJGFhYaHm51QTi8Wwb9++wO8Lq6WlJRw6dKjbYfQclkulSCSCQqHQVLnYDzd/DOzfvz+U5crzpRLLpNKjR486UqmsmdBHR0cxNTWF999/3xsJbmZmxkvwwGYiz2QyeOedd/DGG2/U/DBRFCtq6e7ramPEu7Xy0nXuv6vV6ImIiHarurOtpVIpAKg6ZWqhUEB/fz8SiQQSiUTdD4tGoxVJ2LbtLXutS5IEQRB8ydv9N5M5ERHRpxqenKVaz/JIJAJJknz/1ROPx33PnedyOd9jaZZl+daPjY0hk8l4r+/evYuxsbFGwyYiItoVAs+2tl2pVAqWZcEwDOi6jgMHDvg62WWzWei67r1OJBJwHAeapkHTNPT19TXUGkBERLSb1G1yb4daCbla8z0TOBERUW0dr6ETERFR6zGhExERhQATOhERUQg0ldAXFhbw8ssv4/jx496yO3futCwoIiIiCiZwp7jbt29jdXUVyWTSNzTr6dOnce/ePV+SJyIios4InNBFUUQymQRQOWxrvbHciYiIqD0CN7nXGqGNs6ARERF1R+CEbpomHjx4ULH8wYMHNSdaISIiovYJ3OQ+OjqK8fFxPHjwwBv61bIs9PX14caNG+2IkYiIiOpoaqS4mZkZWJaFxcVFOI4DVVVx+PDhVsdGREREDQqc0BcXF3H48OGGJ2MhIiKi9gt8D/3111/HxsZGO2IhIiKiJgWuoff398MwDBSLRQiCwKZ2IiKiHhA4ob/55pvevwuFAhYWFrBnzx7IsoynnnqqpcERERFRY7Y1lnskEsGJEyfw/PPP47d+67fw+uuvtyouIiIiCiBwDX1tbQ39/f0ANoeBnZ+fx8bGBkZHRxGPx1seIBEREdUXOKGPj48jFoshk8kgHo/jypUrOHjwYDtiIyIiogYFTuiO4yAWi+GNN95oRzxERETUhMD30BOJBE6fPt2OWIiIiKhJe4otnCLtzp07PZfsHz16hHw+jz//8z/3jTX/q7/6qzh37hx+9KMf4atf/WrF+06fPg1VVfHBBx9gbGysYv1Xv/pV/Pqv/zref/99jI+PV6wfGxvD8ePH8Y//+I/4vd/7vYr1Fy9exC/90i8hn8/ja1/7WsX63/3d38Xg4CDu37+Pr3/96xXrv/a1ryEWi+Hv//7vce3atYr1f/RHf4Sf//mfx7179zA3N1exfnR0FIqi4G/+5m/wzW9+s2L93Nwcnn76aei6XnWu+29+85v4mZ/5Gdy6dQt/+7d/W7H+r/7qrwAA3/jGN/Dtb3/bt+7JJ5/EX/7lXwIA/uRP/gS5XM63/vOf/zw0TQMA/OEf/iGWlpZ865955hn82Z/9GQBgamoKKysrvvVf+tKXcPXqVQDAa6+9hu9///u+9QcPHkQqlQIAXLhwAT/84Q+9dYVCAb/8y7+M3//93wew+QP24cOHvvcfOXIEv/M7vwMA+M3f/E18/PHHvvXHjh3Db//2bwMATp06VVE2O+3cM00TP/nJT/CLv/iLAIKde1//+lXs+WwEP/n4Ib7wP30BwOZIk88++2wozr3vfOc7iEQi3vrtnHsAcOjQoR1/7v3cz/0cRFHsye+9bp17oiji1VdfRSwWw759+yre1yp1m9xff/11nD17Fs8//zwA4JVXXqm6XbFYhGmaPZfQiYiIdoO6NfR0Oo1f+ZVf8RL6yy+/jGQy6ftVWrrtzMxMeyJtkltDb/cvo51maWkJhw4d6nYYPYflUunUqVPemBNB/d23/ws+J76AD/77f8ZvnP61NkTXXTxfKrFMKnUqD9WtoSeTSd/rN954Y8sx3Mu3JSIios4I3Cmu1oQsnKyFiIioOwI/traxseEb4tWyLKysrHBcdyIioi4KXEPXdd33WpIknDhxAocPH8a9e/daFhgRERE1bltjuZcrfSyMiIiIOqduk3uhUEAmk0Eul8PGxgYsy4JhGBXbWZYFVVXbEiQRERHVVjehRyIRnDlzBmfOnMH09DSeffZZnD17tmI7SZKqPspGRERE7ReoU9zExAQWFhYqJmMpFAqwbZsJnYiIqEsC30M/ceJExbJIJAJRFNkpjoiIqEsCP7YGbM6Jvri4CMdxfMtXV1dx/PjxlgRGREREjQuc0BcXF6FpGg4ePAjHcSAIAoDNaVVfe+21lgdIRERE9QVO6CsrK3jzzTcBbPZsBz4dIe7BgwfemO9ERETUOYHvoff393v/liQJ77zzjvd6fX29JUERERFRME0PLON2gFteXsbGxgYAVMxLTURERJ3RVC/369ev4+7duwCARCKBF198EV/5yldYQyciIuqSpnq5j46Oev+WJAnvvvsuVlZWKp5PJyIios5o2VjuBw8exI0bN1q1OyIiIgqgZg39lVdeaXhHxWIRpmkGeg8RERG1Rs2EXiwWkUwmGx7SNZ1OtyQoIiIiCqZmQk8mk4HuiyeTyW0HRERERMHVvIcetJPbgwcPGtpO0zRks1noug5d1wN9xtTUVKDtiYiIdoPAvdzv3LlTdbnjONB1ve5Y7tPT0xgYGICiKN7rbDbrva73Xnd0OiIiIvpU4IQ+NzcHWZZ999Udx8Ha2hpUVa37/tu3b2NiYsJ7PTQ0hHQ6XTehm6YZNFQiaqMPPvgAP/zhD/Hxxx9vuY34+QiAf8D+L30RS0tLnQuug8J6XNuxm8rkySefxDPPPIOnn36626EET+iJRAJnzpypum5hYaHme6slZUEQYBhG3c/N5/M4cuQIR6Mj6gEfffQRLMvCl770JTz11FPYs2dPt0Mi6rhisYiNjQ18//vfx5NPPomf/dmf7Wo8gZ9D3yqZA/BmXtuKbdsQRdG3rPx1NdlsFvF4vLEAiajt3n//fTzzzDOIRCJM5rRr7dmzB5FIBM888wz+6Z/+qdvhtG5gGWBznvRayudPb2SdO0VrvR8LRNQ5P/rRj9DX19ftMIh6Ql9fHz766KNuhxG8yb3awDEPHz7ExsZG3XvogiDAtm3fsvLX5TKZTEP35uvJ5/Pb3kfY7Kb7XEGwXPwKhQKAynLZu3dvN8Ih6jl79+7F48ePu/7dETihP3z4EFeuXPEti0Qi6OvrqzsAjSiKFTVx93W1GrhpmpBlOWiIVcViMezbt68l+wqDpaUlHDp0qNth9ByWS6VIJIJCoeArl6WlJTa1E/0r91rY6rvj0aNHHalUBk7oV65c2fL59LW1Nd986eWi0WhF4rZte8ukbdu2r8Pc8vIyLMuCpmlQFAWSJAUNn4iIKJQCJ3Q3mbtzoJfSNA1vvPFGzffH43Hfc+e5XM7XpG5ZFkzThKIokGXZl+x1Xcfa2hoSiUTQsImIiEItcEJfWFjA5OQk+vr6UCwWAWw2NxSLRaytrdVN6KlUCpqmwTAMWJaFAwcO+J5Bz2azMAyj4rl0XdeRzWa9GrqqquwoR0Rt4XbGJdpJAid0y7Lw7rvvVl13/fr1hvZRq4adSCSqrldVtSWd44iI6rl06RKGhoYaGsGylDs0dSqVakdYRDUFTujRaHTLdaOjo9sKhoioFywuLqJQKARO6EG3J2qllj6Hvri42MrdERF1XDabRTKZhGEYNcfOqKa83w9RJwWuoR8+fBh//Md/DGDzYXr3PpPjOLh79y6+9a1vtTZCIqIOsiwLiUQC6XS6ZeNgEHVC4ISeTqdhWRYkScLDhw/x8OFDb507AAUR7U63s3+H/+Pu212N4X8fOoozyotNvbe0M5z7RE61hK5pGqLRKBzHgWEYSCQScBwH6XQaAHDz5k0Am2Np2LYNx3GwvLyMI0eOeDV40zSRTqcRiURw9uxZAJtP/Rw4cIA/IqgpgRP6wMAAksnkluuIiHaqTCbjzRuhKApGRkYqerzrug5JknxN647jIBqNIpFIQNM0b/nk5CTGxsagKAoURcHg4CDu378PYLM/kqqqSKfTkCQJkiQhFothcHCQCZ2aEjih13qU48SJE9sKhoh2tjPKi03XjntBafKWZRmCIFQ0u0uShKmpKW9QrFod4WZmZioGwCr9DPf/7jaltzD52BwFFbhTXF9f35aTsNy4cWPbARERdYNlWd44F+5/sVgMuq77tpNlGalUCrlcDsPDwxgeHt6y85woitA0Dbque6Nels9fwREvqVUC19C/8Y1vYG1tDYVCwRvDHdgc4/3999+vOnkLEVGvMwyj4vlxWZYxPDzs9RtytyvtzT41NQVd16uOnzE8PIyZmZmKx31ZA6d2CJzQC4UCkslkxTzmxWIRc3NzLQuMiKiTVldXK5ZFo1FIkoRsNuslbNM0AcBL6Kqq+uaccLkd4txkXlqLLx0Ns96Mk0SNCpzQk8nklpOznD9/ftsBERF1kmEY0DQN+XweAwMDvnviuq7Dtm2vsuIOOW1ZFrLZLIBPH3MzTRPz8/PI5/PQdR2qqiIej3s94gHg8uXL3tDVpmlC0zRYlgVd1xGPxzE7Owtg82miRCLB5ngKZE/RHZC9Be7du4fjx4+3anct4U5bx+lT/ThNaHUsl0qnTp1CoVDAwsKCt4zlRORX65roVB4KXEO/c+dO1eWO40DX9Z5L6ERERLtB4IQ+NzcHWZYRiUS8ZY7jYG1tjc9OEhERdUnghJ5IJHDmzJmq60qb5IiIiKhzAj+HvlUyB2oPOkNERETt09LZ1rYacIaIiIjaK3CTe7WBYx4+fIiNjQ3eQyciIuqSwAn94cOHuHLlim+ZO2JcaUc5IiIi6pzACf3KlStbDixDRERE3VEzoW9sbFQsYzInIiLqPTU7xeVyOZw8eRKDg4OYnp6uGK/YsiwsLCzg3r17bQ2SiIiIaqtZQz9x4gRs20Y8Hq96f1ySJG+s4Tt37uD06dPtiZKIiIhqqllDv3fvHoaGhhrq7Hb69Okth4UlIiKi9qqZ0G3bxlNPPdXwzlo4zwsREREFUDOhl87f24hCobCtYIiIOi2bzWJ4eBhf/vKXMT097X3v6bqOY8eOeX2IXJZlYWRkBMeOHYOu6zBNEyMjIxgZGan5OVNTU5iammpp7JZlYXp6GoODgxWfb5omLl68WBF/O1iWBU3TMDg4iOHhYWiaVnM50J7y2I5ei6cpxRrS6XSt1RWuX78eaPtO+Pjjj4vf/e53ix9//HG3Q+kp3/3ud7sdQk9iuVR66aWXisePH/ctC1s55fP54nPPPVe0bdu3PJPJFF944YWK7VdXV4u5XM633blz52p+Ri6X872nWCwW5+bmthG1fz8vvPBCcX5+vuq6Tjl69GjVzzt58mTx6tWrvmXVyqObthtPrWuiU3moZg1dFMWGe7AvLi6yyZ2IdqRoNApBEJDJZHzLBUGA4zgVT/iYpglZln3b1SPLsu89ALC8vLyNqP1xXr58GVNTUxUtq27H5W7r6+vzva5WHt3Ua/E0o2ZCHx0dxfz8PN55552aO3nw4AGuX7+O0dHRlgZHRNQp8Xgc2WzWt8xxHMiyXLG8FVrdDK4oCmRZxvj4eEv3SztH3ZHi3njjDYyPj+Pzn/88FEVBLBbz1uXzeWSzWayvr+NP//RP2xknEe0A//W//Sf83dL/2dUYXjz0G/gP//7XAr9PURSMjIzAcRxfjVtRFKTTaaRSKQCbSX6rWq9hGHAcB8vLyzhw4IA3v4Vpmkin0wCAmzdvwjAMWJaFlZUVaJoGQRC8bR3HwezsLAYGBrC8vIwjR440XHOcmZnB0aNHkc1moShK1W0cx4Gu694xWJaFRCLhizMSieDs2bMANscjKT2WVigvDwDQNA3RaNRrEUkkErAsC+l0Gv39/Thy5AhEUawoW3d/tm17ZV9aZo0cU7V4Sv8ObkvNVmXaK+omdEmS8NZbb0HTNMzNzcGyLN86VVVZMyeiHU+WZa/ZXVVVGIbhJYWpqSnvtWEYVb/Y8/m8NzaHLMsYHBz0EkY0GkUikfA6hcmy7DWNu8nUNTw8jLfeeguCIEBRFAwPD+PWrVsNNesLgoBkMonJyUnveMqdO3fOt79sNoupqSmkUilEo1Goqop0Ou0dSywW8x1LPeW3JwD48ka18nB/YJT+cHFbR8bGxpBOp3Ht2jUAmz+wLl68CABeTJOTkxgbG4OiKFAUBYODg7h//773WfWOqTwet5wuX76MaDQKADh27BgkSfJe96KGx3JPJBLeiWdZVs/clyGi3vEf/v2vNVU77hVus7uqqr6autvsXqumLIqi973ovq+8tl9PNpuFIAi+98RiMe9HRiNUVUU2m8WlS5e8JFi6/9L4gM0EOT4+jmQy6fvsZo9FluWKHynlfRPKSZKEqakp2LYNWZYrfjCVDzl+9uxZjI+Pe2UyMzNTkZNK4w16TG4LSmnynpmZ6elkDjQ5HzqTORGFkaIoXrN5+fJMJlMzqbXie9GtyRqG4f3n3hsPIpVKYWFhoeLev2VZEEWxYntBECpaX7fiPrLn/letRh6ULMtIpVLI5XIYHh7G8PBwzcemJUnyrRdFEZqmQdd1Lx7btive06hqldZeT+ZAE7OtERGFldtMfenSJd800fF4HFNTU7576a3kJlM3UW23t7UkSV7TezKZ9C0vT3RA7X4B5dx7zK3k3s4ovcWh63pFTd9VnnCHh4er1qCDtpC4JEmquE2wEzRVQyciCqvDhw+jUCj4EoEgCHWTbLVEWYskSd5gXG6CUhSlorZsWVbdWvDq6mrFskQiUVGTVRQFjuP49p/NZnHixAnf8QY9lkasr69vuc40Td8xlt9eWFlZ8b3WNM3Xkc+2bS+Zlx5v6T6DHJMsy5Akyfd+0zRhmmbD++gG1tCJiEqcPXu2au1MUZSqtVjTNKFpGizLgq7riMfjmJ2dBQCk02kkEgk4joP5+Xnk83noug5VVRGNRiFJEnRd9+3v1q1bXu9qoPaPCXekuMXFRQDAxMSEb/3MzExFs/tbb73l7d/tGe7ea2/kWKqVgWVZyGazsCwLmUzG67Vfuhz4NBGbpukrD/dHjBtrac97AF45iaIIy7Kgqqp3nz0ajSIej3u95AHg8uXL0DQNqqo2/fe5desW0um0d5uikR913banGPLRYB49eoR8Po9YLIZ9+/Z1O5yesbS0hEOHDnU7jJ7Dcql06tQpFAoFLCwseMtYTtQp2WwWd+/erejg12tqXROdykNsciciop7GeUIaw4ROREQ9yTAM6LruNYVTbbyHTkREPSkM46t3EmvoREREIcCETkREFAJM6ERERCHAhE5ERBQCXekUp2mabwjCWpMOuFP9AfCmxWvlNH5ERERh0PGEPj09jYGBAW+Un+np6Zpz987OzvpGPzp27BiA2j8CiIiIdpuON7nfvn3bl7yHhoa2fL6wfMxhYDORl85ZS0RERB1O6NUGthcEoebEA4uLi76kXj5xAREREXW4yd227Yq5eKvNzesSBAH379/3LcvlchxogIiIqExHE3qtCesbmbfWcRwsLi7i1q1bgT87n88Hfk/YLS0tdTuEnsRy8XPH0Wa5ENXW7WukowldEISKOWmDzFF76dKlqpPYN4KzrflxtqzqWC6VIpEICoWCr1y6/cXVStlsFnNzczBNE6Ojozh//jwEQYCu69A0DbZt48yZM17nXMuyMDU15U3xGYvFkE6nAQA3b97c8nOmpqYAAKlUqmWxu1OOzs3NQZIkxONxJBKJLZe3K47t6LV4tqPebGvt1tGELopiRS3dfV2vdq5pGs6ePcvmdiJqKXee8+HhYS+ZA5sdcEVRxOTkpO9JG0mSkEqlYFmW932kqmrdyUOqPcnjzg/eLEmSkEgkvHm+3X25yzOZDA4fPuz7jK2eKOqWXotnJ+top7hoNFqRuG3brpuks9ksotGot12tTnREREG5302ZTMa3XBAEOI5T8Z1jmqbve6tehQSoPtHI8vLyNqJuTF9fX904uqnX4tnJOv7YWjweRzab9V7ncjnfM+VuU5HLMAzYto1YLOY9xlattzwR0XaUfzcBmy2IsixXLG+F6enplu+TdreODyyTSqWgaRoMw4BlWThw4ICvySWbzcIwDCiKAsdxMDIyAuDT+ywAcOLEiU6HTUQN+OD//n/wwf+10tUYnv7fDuLpf/dvA79PURSMjIxUdNBVFAXpdNq7x+s4DiRJqroPwzDgOA6Wl5dx4MABr7JimqbvPrv7/beysgJN0yAIgret4ziYnZ3FwMCANzpmq2qw5XEAm83+0WjUa4lw78Gn02n09/fjyJEjEEWx4pjc/dm27R1zaazuZ0UiEZw9exbAZgWuVrmUH7/bQsJm+cZ0ZejXWveMEomEt14QBPzDP/xDp8Iiol1MlmWv2V1VVRiG4SWnqakp77Vb4SiXz+chSRIkSYIsyxgcHPQSVzQaRSKR8AbFkmXZ6z9U/n04PDyMt956C4IgQFEUDA8P49atW3Wb9avdiiwfs6M8Dl3XvXhdbqvE2NgY0uk0rl27BmDzh83FixcBfDpS5+TkJMbGxqAoChRFweDgoPeocTQahaqqSKfTXrnEYrGa5QIA586dw+XLl73Oz8eOHYMkSU11ht5tupLQiSicnv53/7ap2nGvcJvdVVX11dTdZvdaNWVRFL2au/u+Rh7HLZXNZiEIgu89sVjM+5FRiyzLFT8OyvsElJMkCVNTU15fpvIfKgcPHvS9Pnv2LMbHx71YZmZmKlorSo/Z/X+j5eK2XJQm72afbNqNONsaEdG/UhTFazYvX57JZGom6K2a4YNwa9SGYXj/KYrStk5jsiwjlUohl8theHgYw8PDNccLkSTJt14URWiaBl3XvRaC8keRg5SLZVkV2zOZN44JnYjoX7nN7pcuXfIl0Xg8DsdxkE6n25JcLcvykpnb5F36Xyt+LFTj3ka4du0a7t+/j1gsVvPxu/KEOzw8DFmWoapqRbN9MyRJ4tDe28CETkRU4vDhwygUCr6auCAIdRN5kEGygM3k5Y7C5yZKRVEq5quwLGtbj+qur69vuc40Td++y5v1V1b8HRxLn5t3O8S5NejSJF66zyDl4v54KX2/aZp8sqlBvIdORFTi7NmzVWuJ7gA05UzThKZpsCzLG+BldnYWAJBOp5FIJOA4Dubn55HP56HrOlRVRTQahSRJFTXiW7dueb28gdo/JtzHfC3LQiaT8XrLly4HPk3Epmn64nB/PLiP5bmj37nc+ERRhGVZUFXVu88ejUYRj8e9XvIAcPnyZWiaBlVVmy6XW7duIZ1Ow7IsiKLY0I8p2rSnWCwWux1EO7lD7nHoVz8OcVody6XSqVOnUCgUsLCw4C1jOYVfNpvF3bt3vV7uVFuta6JTeYhN7kREVJV7S4B2BiZ0IiLyMQwDuq57TeG0M/AeOhER+XB89Z2JNXQiIqIQYEInIiIKASZ0IiKiEGBCJyIiCgEmdCIiohBgQiciIgoBJnQiIqIQYEInIiIKASZ0IiKiEGBCJyIiCgEO/UpEu5plWZifn8f169e9KUGBzXnELcvC2bNnvWFQTdNEOp0GANy8eXNbn2sYBtLpNPr7+zmjGbUEEzoR7WqSJGFiYgILCwtQVRWqqnrrHMfB4OAgbt68CVmWEY1GkUgkoGnatj9XlmWoqopcLrftfREBbHInIgIACIJQdZksy5ifn2/LZ4qi2Jb90u7EhE5EVEM+n8fAwEC3wyCqi03uRNQy//w/HuP//R+PuxrD//yFvfg3X9i77f1YlgVN0zA2NoZEIrHldqZpwrZtOI6D5eVlHDlyxDf1qOM4mJ2dxcDAAARBgOM4UBSlYj+O4+Do0aOIxWJIJBKcvpQCY0InIvpXpmkim816ryVJqtoUX2pychJjY2NQFAWKomBwcBD379/31p87dw6XL19GNBoFABw7dgySJHmvXbZtI5lM+u7hEwXBhE5ELfNvWlQ77pZoNFpRex4ZGYFpmkilUlXfMzMzA0mSfMscx4EgCDAMA5Zl+ZL3zMxMRTI3TROGYdRsCSCqh/fQiYhqUFUVuq5vuV4URWiaBl3XYRgGgM3aNrDZbF+e7MuT+crKCvL5PObm5uA4Toujp92ENXQiohrc5LyV4eHhqrVux3EgSRIsy6r5fvfxNcdxcOnSJT6TTk1jDZ2ICNiydqzr+pb3td0OcW4yL92HYRiQZRmSJHk1d/c9pmlW7CuRSGBlZcV3D58oCNbQiWhXsywL2WzW+7+blNfX17GysoLDhw9jYmICwGYynp+fRz6f9xJ9PB6HpmleUr98+TI0TfN+BNy6dQvpdBqWZUEURe/ZdsMwoOu697mKokAQBIyPjyOZTPJ+OgW2p1gsFrsdRDs9evQI+XwesVgM+/bt63Y4PWNpaQmHDh3qdhg9h+VS6dSpUygUClhYWPCWsZyI/GpdE53KQ2xyJyIiCgEmdCIiohBgQiciIgoBJnQiIqIQYEInIiIKASZ0ImrKJ5980u0QiHpCr1wLTOhEFNgTTzyBjz76qNthEPWEjz76CE888US3w2BCJ6Lgnn32WXzve9/DxsZGz9ROiDrtk08+wcbGBr73ve/h2Wef7XY4HCmOiIJ7+umnAQA/+MEP8OMf/7jL0RB1zxNPPAFJkrxropuY0ImoKU8//XTdL7G/+/Z/wefEF/DBf//P+I3Tv9ahyDqHI+ZVYpl0D5vciYiIQqArNXRN0yBJkjct4VYzGTW7PRER0W7T8Rr69PQ0JEmCoihQVRWrq6s1pwsMuj0REdFu1PGEfvv2bSiK4r0eGhqCrust256Iumvp3fv4T3+Zwd3/+G381PnZbodDtGt0tMndNM2KZYIgwDCMlmxPRN23+o//H37uy7+E4ic/BX76E/z0Ixtf/F++2O2wiEKvowndtm2IouhbVv56O9tX4073zkdrKj169KjbIfQkloufKIr47Gc/23C5fHbvZ7EHjxH7X4vYv/8pAE8B+EJoyzWsx7UdLBM/N/+4+ahdOprQHcepuU4QhG1tX83jx48BAO+9916DUe4e+Xy+2yH0JJaL36uvvgqg8XJ59ktfAPB9/OAHbQyqh/B8qcQyqe7x48d48skn27b/jiZ0QRC8nuqu8tfb2b6a/fv347nnnsPevXuxZ8+eQO8lIiLarmKxiMePH2P//v1t/ZyOJnRRFCtq3e7rarXtoNtX85nPfAaRSKSZcImIiFqinTVzV0d7uUej0YpEbNs2ZFluyfZERES7VccfW4vH477nyHO5nG+gGMuyfOvrbU9ERETAnmK7u91VoWkaotEoLMsC4B/5TdM0GIaBmzdvNrQ9ERERdSmhExERUWtxchYiIqIQYEInIiIKASZ0IiKiEOjK9KnUHY7jQNd1rK+vY2Jiou72u2Xa2iDHqes6TNP0JgzKZrNIJBKQJKkjsbYLpzSujueGH79Dehtr6LuEYRgwDAOWZaFQKNTdfrdMW9vMcWYyGYyMjCCdTkNV1R3/hc0pjavjueHH75AdoEi7ytWrV4uTk5N1t3vhhRd8r/P5fPHcuXPtCqtrgh7n/Px8u0PquKBlwHOjujCeG9XwO6R3sYZOFXbLtLW75Thr4ZTG1e2W42wXll938B46VWjFtLU7QbPHqes6RFEMxX3BbkxpvBPw3Nie3XKe9BomdKrQimlrd4JmjjMWi0EQBO/e6MWLFyGKotcRaqfpxpTGOwHPje3ZLedJr2GTO1VoxbS1O0EzxxmNRn0dnQYGBjA3N9eW+DqhG1Ma7wQ8N7Znt5wnvYY19B0qm83i7t27Nbfp6+tDKpUKvO9WTFvbLUHKpZnjNAzDN9ufJElV7xfuFN2Y0ngn4LmxPbvlPOk1TOg7lKIobWvK28nT1gYpl6DHaVkWRkZGcP/+fd/7dvKjSZzSuDqeG9uzW86TXsMmdwKwe6etDTKdryRJSCaTvi+qu3fvIpFIdC7gNuCUxtXx3Ahmt54nvYSzre0SpmnCMAzoug5gs/etLMuIRqMAdve0tUGm8y390lpfX8eBAwdCUS6c0rg6nhuf4ndI72NCJyIiCgE2uRMREYUAEzoREVEIMKETERGFABM6ERFRCDChExERhQAT+hZqjUVMtNOE9XwO63ERNXNu9/xIcZqm+Z57VBSl7aMv6bruPS+paRrm5uYgiqK3bH19HYVCAdFotOK51K22jUQimJiYCBRH6XOtpc97dlu14wQ241VVtSNx6roOTdOQTCa7NvmFe24KgoB4PO4ttywLhmFAEAS89dZbXYmtXCaTQTweb9uwm7v5OgV681otPT/Lz0P3+nEHxGlFvKZpIp1OA4DvWfSdohvncC1NXbPdnY69MRcuXChevXq1I581Pz9fXF1dbejzJycni5OTkw1tm8lkiufOnQsUy8mTJ4u2bRfn5uaKc3Nzgd7bblsd58mTJ4uZTKYjMUxOTnbss7ayVTnYth34790qtm1XXd7ua2i3XqfFYu9eq/Pz88WjR48W5+fnq67b6lxpVj6fryi/q1evVpR/r+rEOVytPFp1zbLJvYTjOFhdXa36q6yvr69iWTKZhK7r3ihItbZ1xxifmppqOBZ3msFEItGTQ0hWO86xsTFMTk525PMjkUhHPqcZgiB0bdzqTCZTdfmRI0d8Q3HuVL10nbrx9Oq1KooiUqmUV3MuJUlSy1tsqu1vaGiII8SVqFYerbpmmdBLZDIZDA0NNby9e/I2OqOSqqrIZDIVXyzV2La9Iyd2EATB+4Lb7aLRaFfKIZfLVV0uy3Ldmeh2gl66ToHev1ZlWUYsFsP09HRXPj8ajfbELYheUa08WnXNMqGXyOVygU4898s6yHtisVgoaklbyeVykGV5V0+R6CaOWCzW8c+enp5GoVDYcn1fX1/DiapX8ToNLpVK4fr16zv+bx9Grbxme75TXBC6rkMURQCbv5rLmzXc9ZZleRd3NptFKpVqqiY1OzuLEydOBPp1LkkSlpeXa25jGAZyuRwsy/I6rkiShMnJScRiMSiKAsuykMvlcO3aNd+xAf5jdzuqRCIRnD9/HrZtw3EcLC8vY2JiAoZheNtJktR0BzPHcaDrOgqFQtUOMe6Xo1v2bnN0vfhKucfotgJUuwi6XQ6WZcG2bQD+5sdacdX6u7qTWziO43ufruveeeeuE0URhULBO2+AzdpmaRyyLMMwjK42gYblOgV2zrUqSRJUVcXU1FTNDmvNxrYVy7K82xc3b94MtJ+tzn1Xre+UWuXfCs2ew4lEwlce2Wy2tdds03f2O6iRjgoXLlzwdZIp75Q0OTlZzOVy3rqjR4/63p/P57f8jAsXLhQvXLhQzGQyxUwmU5yfny9evXq1akeTCxcu1OwUc/Xq1eLJkydrHkuxWCyurq5WdC7J5XLFkydPFldXV70OOO5n1jr20veVxll+vC+88ELduErff+HChWIulytmMpni3Nxc8erVq1U7d5SXx7lz53yxbBWf+/cqFotVy7u8A163yuHcuXPe8R89etQXd9C4yv+u+XzeVwbu+Ve+P/cYqnVKKlXrPN+u3XidFou9fa2WXh+2bRdfeOEFr3y3c56Wvqd0P9XKovycbGQ/W537riDfKaXlX0+nzuHS7Vt5zYaiyd00Taytrfl+gQuCAEmSvMcQMpmM1wRa7Z6abdtVO8m43F/EiqJAVVVMTEw0Xcvp7+9v6n2iKMJxHK8zSyKRaOjYS99Xejxb7b9RAwMDkGUZiqIgkUhgaGgIw8PDFdstLy/7mi8PHjzo1TZqxec2MzmOg9u3b1eUd2mTdjfL4eDBg0gkElXPiaBxuX9Xy7KwsrLiayYeGhry3jM/P+/bX6P3lAVB6Fqz6265ToHeu1bdz0omk1U7/G0ntnrnU/ntt3r7qXfuA8G+U9zyb4VWnMNBb0cGuWZ3fJO7YRiwLKvqxSdJkleQkiTBtm2vMMs7sliW1ZH7vpZlbasDTfl78/l83WMHqp9Etb4YmxGNRmHbNgzD8PXwLm3qsiwLhUKh4suoVtkbhlG3zHqlHFRVRT6fDxxX+fG5z7GXfkk5juM963/x4kV8+ctf9n5QNZq03Gb5Tttt1ynQm9eqqqre89al8W0ntmbUu963OvddjXyntLqjYqvO4aCCXLM7PqHX61G9vr4OYPNEnp+fx/nz55HJZDA2NuY7qURR9O59ttPi4iLOnz/f9PvLH9Vq5NgBePd7SrXji9G9b1TKvT/k3ueq9rhZtfiC6JVyKH9crdG4qv1d+/v7Kx59UxQFjuPg2rVrcBwH+XwemqbBNE2kUqmKzyhPTI1+uYyMjPiuh3g8vq1azm67ToHevVZTqRTGx8eRTCZbElszau2n1rnvauQ7pdWPtbbqHK6n2WsW2OEJ3f3VJsty1ef4LMvCkSNHAGyeQOfPn0c+n686+o4gCA11gtmObDaLw4cPt/QRjkaOvZMEQcDq6qpv2fDwMN5+++2qJ3UjzYbRaLRuk1OvlYOr2bii0eiWz6bOzs5iYmLC+/EgyzJGRkaqbut2oHI5jtPQF10rR/ridbqpV85R9zE2TdO8lp1eiQ2ofe67tvudElQrz+F6mr1mgR3+2NqlS5cgCAKi0Sj6+/t9TUNu7cU9Yd3mm60eqYrFYi27t1j6i9ZlmiZ0XceVK1e2te/yppdGjh1A1VrNdk/8ascZi8WwsrICYPOY3377bQD+GoZ7DG7v1XrxSZKEeDzua4IDNv+m7vu6WQ61NBpX+d9VlmWIoljx7LR7n670fiLwafNi+T3N8nPdsiwMDAxs44iC243XKdAb1+pWP35SqZSvHDsdW6391Dv33eX1vlNaeWupledwuVZesz1fQ9c0DYuLi1hZWfHuI7ljZVuW5XUGunbtGjRN83WsuHXrlrcfRVEwPDzsPfIEACdOnPAelRAEoeoJoGkaVlZWvM+v1ezoxrq2tuY9glCq0VqPZVmYnp72mlNVVYVlWZidnfUtc//wtY7dNE3Mzs56TVSJRALZbBYLCwveCSrLsvf+dDqNRCJRs4nHHaCidJ+uZDKJS5cueR1WFEXBmTNnvOYxURSRTCaRTqe9e+ONxJdKpaBpmjcil+M4kGXZuxcoy3LHy6H83Nhq7OdG4qr2d7158yY0TUM+n/eaKN17oJIkeV8cjuN4fwNBEKCqqvfYTPmjTYZh+JpaW2U3Xqdu/L16rU5NTXmPkiaTSV+ikCQJo6Ojvu23G5skSb6ycDt3ust0XUcsFmvoGLc694HNHx+NfKdUK/9aOnEOl5eHG1urrtk9xWKx2NCWO5hpmjAMo+Ii1zQN6+vr3peFpmk9M7ECUTtcvHixpc/jthKvU9rpGj2Hgwhyze7oJvdG6bpedRCGRCLhNQ8Dm78AwzA0JlE1hmEEGjK103id0k7X6DncqKDX7K5I6O5IO+XKH68SBAEHDhzg8IgUSrlcrmtTzTaC1yntdI2ew40Kes3uiiZ34NNnCEuHywRQtbDK7wsT7XS6rrd1PvRW4XVKO12Qc7iWZq7ZXZPQg3I7XxGFQVjP57AeF1Ez5zYTOhERUQjsinvoREREYceETkREFAJM6ERERCHAhE5ERBQCTOhEREQhwIROREQUAv8/oc4CWRlzPRQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 540x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# utility for making cdf plots (diff of level or log PD)\n",
    "def cdf_pd_diff(df,race,plotrace,clfs,lim = (-1,1),log=False):\n",
    "    fig, ax = plt.subplots(1,1,figsize=(7.5,5))\n",
    "    if log:\n",
    "        diff = df[clfs[1]].apply(np.log) - df[clfs[0]].apply(np.log)\n",
    "        factor = 1\n",
    "    else:\n",
    "        diff = df[clfs[1]] - df[clfs[0]]\n",
    "        factor = 100\n",
    "    for group in plotrace:\n",
    "        x = np.sort(factor*(diff[race==group]))\n",
    "        y = np.linspace(1,x.shape[0],x.shape[0]) / x.shape[0]\n",
    "        ax.plot(x,y,label=group)\n",
    "    ax.set_xlim(lim); ax.set_ylim((0,1))\n",
    "    ax.set_xticks([-1,-0.5,0,0.5,1])\n",
    "#     ax.set_yticks([0,0.25,0.5,0.75,1])\n",
    "    ax.axvline(0,color='k')\n",
    "    ax.axhline(0.5,color='k',linestyle='--')\n",
    "    plt.xticks(fontsize=16); plt.yticks(fontsize = 16)\n",
    "    if log: \n",
    "        ax.set_xlabel('Log(PD from %s) - Log(PD from %s)' %(outnames[clfs[1]],outnames[clfs[0]]),fontsize = 18)\n",
    "    else:\n",
    "        ax.set_xlabel('PD from %s - PD from %s' %(outnames[clfs[1]],outnames[clfs[0]]),fontsize = 18)\n",
    "    ax.set_ylabel('Cumulative Share',fontsize=18)\n",
    "    ax.legend(frameon=True,framealpha=1,fontsize = 16,loc='lower right')\n",
    "\n",
    "# make figures\n",
    "race = df0['Race']\n",
    "cdf_pd_diff(df0,race,plotrace,clfs,lim = (-1,1),log=False)\n",
    "cdf_pd_diff(df0,race,plotrace,clfs,lim = (-1.5,1.5),log=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Table VI: Decomposition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>roc</th>\n",
       "      <th>precision</th>\n",
       "      <th>brier_score</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>model</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Logit</th>\n",
       "      <td>0.498538</td>\n",
       "      <td>0.007465</td>\n",
       "      <td>0.007415</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>LogitNonLinear</th>\n",
       "      <td>0.501351</td>\n",
       "      <td>0.007496</td>\n",
       "      <td>0.007415</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>RandomForest</th>\n",
       "      <td>0.498721</td>\n",
       "      <td>0.007448</td>\n",
       "      <td>0.007451</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>RandomForestIsotonic</th>\n",
       "      <td>0.499342</td>\n",
       "      <td>0.007448</td>\n",
       "      <td>0.007415</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                           roc  precision  brier_score\n",
       "model                                                 \n",
       "Logit                 0.498538   0.007465     0.007415\n",
       "LogitNonLinear        0.501351   0.007496     0.007415\n",
       "RandomForest          0.498721   0.007448     0.007451\n",
       "RandomForestIsotonic  0.499342   0.007448     0.007415"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>roc</th>\n",
       "      <th>precision</th>\n",
       "      <th>brier_score</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>model</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Logit</th>\n",
       "      <td>0.500043</td>\n",
       "      <td>0.007489</td>\n",
       "      <td>0.007415</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>LogitNonLinear</th>\n",
       "      <td>0.501351</td>\n",
       "      <td>0.007496</td>\n",
       "      <td>0.007415</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>RandomForest</th>\n",
       "      <td>0.497118</td>\n",
       "      <td>0.007426</td>\n",
       "      <td>0.007454</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>RandomForestIsotonic</th>\n",
       "      <td>0.499291</td>\n",
       "      <td>0.007463</td>\n",
       "      <td>0.007415</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                           roc  precision  brier_score\n",
       "model                                                 \n",
       "Logit                 0.500043   0.007489     0.007415\n",
       "LogitNonLinear        0.501351   0.007496     0.007415\n",
       "RandomForest          0.497118   0.007426     0.007454\n",
       "RandomForestIsotonic  0.499291   0.007463     0.007415"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "adding race first\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Race</th>\n",
       "      <th>Technology</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>roc</th>\n",
       "      <td>0.0</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>precision</th>\n",
       "      <td>0.0</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>brier_score</th>\n",
       "      <td>0.0</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             Race  Technology\n",
       "roc           0.0       100.0\n",
       "precision     0.0       100.0\n",
       "brier_score   0.0       100.0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "adding technology first\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Technology</th>\n",
       "      <th>Race</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>roc</th>\n",
       "      <td>97.52</td>\n",
       "      <td>2.48</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>precision</th>\n",
       "      <td>146.89</td>\n",
       "      <td>-46.89</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>brier_score</th>\n",
       "      <td>100.37</td>\n",
       "      <td>-0.37</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             Technology   Race\n",
       "roc               97.52   2.48\n",
       "precision        146.89 -46.89\n",
       "brier_score      100.37  -0.37"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# import scores\n",
    "score00 = pd.read_csv(path + '../output/eval_output_race0_interestrate0.csv',index_col=0)\n",
    "score10 = pd.read_csv(path + '../output/eval_output_race1_interestrate0.csv',index_col=0)\n",
    "display(score00)\n",
    "display(score10)\n",
    "\n",
    "# decompositions using saved scores\n",
    "x_part = score00.loc[clfs[0]] - score00.loc[clfs[1]]\n",
    "total = score00.loc[clfs[0]] - score10.loc[clfs[1]]\n",
    "g_part = score00.loc[clfs[0]] - score10.loc[clfs[0]]\n",
    "\n",
    "racefirst = (100*pd.DataFrame({\"Race\": g_part / total, \"Technology\": 1 - g_part/ total})).astype('float').round(2)\n",
    "\n",
    "techfirst = (100*pd.DataFrame({\"Technology\": x_part/ total , \"Race\": 1 - x_part / total})).astype('float').round(2)\n",
    "techfirst = techfirst[[\"Technology\",\"Race\"]]\n",
    "\n",
    "print('adding race first')\n",
    "display(racefirst)\n",
    "print('adding technology first')\n",
    "display(techfirst)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
